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Abstract 

We perform a detailed study of double Higgs production via gluon fusion in the Effective Field 
Theory (EFT) framework where effects from new physics are parametrized by local operators. Our 
analysis provides a perspective broader than the one followed in most of the previous analyses, 
where this process was merely considered as a way to extract the Higgs trilinear coupling. We 
focus on the hh 6677 channel and perform a thorough simulation of signal and background 
at the 14TeV LHC and a future 100 TeV proton-proton collider. We make use of invariant mass 
distributions to enhance the sensitivity on the EFT coefficients and give a first assessment of the 
impact of jet substructure techniques on the results. The range of validity of the EFT description 
is estimated, as required to consistently exploit the high-energy range of distributions, pointing out 
the potential relevance of dimension -8 operators. Our analysis contains a few important improve¬ 
ments over previous studies and identifies some inaccuracies there appearing in connection with the 
estimate of signal and background rates. The estimated precision on the Higgs trilinear coupling 
that follows from our results is less optimistic than previously claimed in the literature. We find 
that a ~ 30% accuracy can be reached on the trilinear coupling at a future 100 TeV collider with 
Sab^ 1 . Only an 0(1) determination seems instead possible at the LHC with the same amount of 
integrated luminosity. 



I. INTRODUCTION 


Unveiling the dynamics at the origin of clectroweak symmetry breaking (EWSB) is a 
primary goal of the experiments performed at the Large Hadron Collider (LHC) at CERN. 
The discovery at Rnnl of a new boson with mass m/, ~ 125 GeV and properties similar to 
those predicted for the Standard Model (SM) Higgs boson has been a major leap forward 
in this direction DEI- The many direct searches and precision measurements performed 
by the ATLAS and CMS collaborations all indicate the existence of a gap between the 
clectroweak (EW) scale and the scale of New Physics (NP), unless the latter is very weakly 
coupled to the SM sector. This justifies the use of Effective Field Theory (EFT) to give a 
low-energy parametrization of NP effects in terms of a series of local operators. Measuring 
the coefficients of these effective operators would give access to a wealth of information 
on the UV dynamics, most importantly whether it is strongly or weakly coupled. Higgs 
studies at Runl have mostly focused on on-shell single-production and decay processes, thus 
probing the strength of the underlying EWSB dynamics at the scale Q = nth■ They set 
limits on possible modifications of the Higgs couplings, whose naive expected size is of order 
8c/c ~ {gl/g‘sM) m h/ rn l^ where m* is the mass of the new states, g * is their coupling strength 
to the Higgs boson and gsM a SM coupling. The improved performances in energy and 
luminosity which will characterize the LHC Run2, on the other hand, give the opportunity 
to directly probe the EWSB dynamics at much higher energies (Q ~ E nth) through the 
study of 2 —)■ 2 scattering processes. For a typical scattering energy E , effects from New 
Physics are expected to be of order ~ (<?*/ g 2 s M )E 2 /m%, hence enhanced by a factor E 2 /m\ 
compared to those entering on-shell Higgs processes. Exploring higher energies thus gives 
access to potentially larger corrections, but at the same time poses the issue of assessing the 
validity of the EFT description. Determining at which point this latter breaks down in fact 
requires adopting a power counting to estimate the size of the local operators in terms of 
the parameters (masses and couplings) characterizing the UV dynamics. At the same time, 
the power counting puts the limits on the effective operators into perspective, and helps 
inferring how much theoretical space is being probed. 

Double Higgs production via gluon fusion is one example of scattering process which can 
disclose key information on the EWSB dynamics, its underlying symmetries and strength. 
It is the only process potentially observable at the LHC that can give access to the quartic 
couplings among two Higgs bosons and a pair of gluons or of top quarks, as well as to 
the Higgs trilinear self-coupling. Previous studies of this reaction mostly focused on the 
extraction of the trilinear coupling in the context of the SM mm- An analysis based on the 
wider EFT perspective, however, can give access to a much richer spectrum of information 
on the UV dynamics. In this work we provide such analysis at the 14TeV LHC and a 
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future 100 TeV proton-proton collider by focusing on the hh — » 6677 final state. We perform 
a detailed Montecarlo (MC) simulation of signal and background, and use the kinematic 
distributions to maximize the sensitivity on different effective operators. Previous studies 
in the EFT context appeared in Refs. [121 IS], and a first analysis of the impact of New 
Physics on the kinematic distributions can be found in Ref. m 

The paper is organized as follows. Section [XT] contains a detailed discussion of the 
parametrization of double Higgs production within the framework of EFT. The power count¬ 
ing of the relevant coefficients is reviewed and a few scenarios are examined where NP gives 
large modifications to the Higgs trilinear coupling whereas other couplings stay close to their 
SM values. The relevance of dimension -8 operators is also investigated and the validity of 


the EFT description is assessed. Section III reviews the phenomenology of double Higgs pro¬ 
duction at the LHC and a future 100 TeV collider, giving a first estimate of the sensitivity 
on the kinematic tail at large invariant masses of the Higgs pair. Our Montecarlo analysis 
of the 6677 decay mode is illustrated in Section IV, while results on the sensitivity to the 
coefficients of the effective Lagrangian are collected in Section [Vj We compare with previous 
analyses of the 6677 final state in Section VI Section VII puts our study into context by 
briefly discussing existing studies of other final states and gives an outlook on possible future 
developments. We collect useful formulas and further details on our simulation of the 6677 
background into the Appendices. 


II. EFFECTIVE PARAMETRIZATION AND POWER COUNTING 

Corrections due to the exchange of new heavy states can be conveniently parametrized 
by means of low-energy effective Lagrangians. There are two formulations which are suited 
to the study of Higgs physics. The first one assumes that the Higgs boson is part of a 
weak doublet, as in the SM, and that SU(2)l x f/(l)y is linearly realized at high energies. 
It thus referred to as the “linear” Lagrangian. In the second, more general formulation, 
SU(2)l X 1 /( 1 )y is non-linearly realized, hence the name of “non-linear” Lagrangian, and 
the physical Higgs boson is a singlet of the custodial symmetry, not necessarily part of a weak 
doublet. Both parametrizations have been reviewed in Ref. [T3j, to which we refer the reader 
for a more in-depth description. The experimental data collected at the LHC during Runl 
seem to indicate that the couplings of the newly discovered boson have values close to those 
predicted for the SM Higgs. Although such experimental information is still preliminary and 
awaits the confirmation of Run2, it motivates the use of the linear Lagrangian for future 
studies. Indeed, small deviations from the SM can be obtained if the Higgs boson belongs 
to a doublet, provided the new states are much heavier than the weak scale. The non¬ 
linear formulation is still useful, however, in those cases where large deviations in the Higgs 
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couplings are allowed. As we will see later in our analysis, this is especially true for double 
Higgs production, where additional couplings not accessible via single Higgs processes can 
be extracted. 

In the linear Lagrangian, the operators can be organized according to their dimension: 


Crlin — £sM + + A£§ + ... 


( 1 ) 


The lowest-order terms coincide with the usual SM Lagrangian C S m , whereas C n contains 
the deformations due to operators of dimension n, with n > 4. 1 For our purposes it is 
sufficient to focus on the operators involving the Higgs boson. The relevant ones in C§ are 2 


A£ 6 D 


^d„(H'H)d“(H'H) + (H'Hq L HH R + h.c.) 
v z zir v ' m^ v ^ 


( 2 ) 


where v = l/(\/2 Gp) 1 ^ 2 = 246 GeV and = 125 GeV is the physical Higgs mass. In order 
to classify the various operators in our effective Lagrangian and estimate their coefficients we 
adopt the SILH power counting [16]. This is based on the assumption that the NP dynamics 
can be broadly characterized by one mass scale m*, at which new states appear, and by one 
coupling strength g*. The latter in particular describes the interactions between the Higgs 
boson and the new states. When building higher-order operators starting from the SM 
Lagrangian, each extra insertion of the Higgs doublet is weighted by a factor 1// = g*/m *, 
while each additional covariant derivative is suppressed by m*. If the Higgs is a composite 
Nambu-Goldstone (NG) boson [T7], the operators Oq, O u and O g can be generated only 
through some small explicit breaking of the global invariance since they violate the Higgs 
shift symmetry [[16]. 3 In this case the SILH estimate of the coefficients in Eq.Q is: 


CH ? Cu i ^6 rN_/ 




(3) 


where A is a weak spurion coupling suppressing c 9 , while the suppression of c§ and c u is 
controlled respectively by (m 2 / 2v 2 ) and yt . In realistic models, the minimum amount of 
explicit breaking entering c g is given by the top Yukawa coupling y t , so that one expects 


1 Here and in the following we assume that baryon and lepton numbers are conserved by the NP dynamics 
and thus by the series of higher-dimensional operators. The only dimension-5 operator invariant under 
the SM gauge symmetry violates lepton number and can be omitted. 

2 For simplicity we focus on CP-conserving operators, CP-violating ones can be included in a straightforward 
way. We omit the operator \H^f)^H\ 2 since it violates the custodial symmetry and is strongly constrained 
by LEP data. Its inclusion has no impact on our analysis. 

3 We denote as O,; the operator in the linear Lagrangian whose coefficient is Cj. 
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Dt A < g*. For example, one has A — y t in models with fully composite £r, whereas partial 
compositeness of both ti and tn leads to A = yjg*y t [16] . 

In the case of the non-linear Lagrangian, the terms relevant for our analysis are 


non—lin —) tt 


c h + c h ' 2 } c<h 3 + 9 * (c h +c ^ W G a ^ 
% + J - <*— h + ^ ^c g - + c 2g — j G^G , 


(4) 


where h denotes the physical Higgs held (defined to have vanishing vacuum expectation 
value). Compared to the linear Lagrangian, the couplings c t of Eq. Q effectively resum 
all the corrections of order {y 2 /p). Therefore, the non-linear Lagrangian only relies on the 
derivative expansion, in which higher-order terms are suppressed by additional powers of 
(E 2 /ml). The linear formulation further requires (v/f) < 1 for the expansion in powers of 
the Higgs doublet to be under control. When both parametrizations are valid, the coefficients 
of the two Lagrangians are related by the following simple formulas: 


_ i °h - 

Cfc 1 2 ^ u ? 


1 , _ . 3_ 

C2i — — 2 V ° H 4“ > c 3 — 1 — pH + Cq , 


= C25 = (S) ’ (5) 


where a 2 = g 1 / 47T. Notice in particular that the same operator O u which gives a (non- 
universal) modification of the top Yukawa coupling also generates the new quartic interac¬ 
tion tthh. 

It is worth at this point to comment about the normalization of the operator Oq in Eq. (|2]). 
This differs from the one of Ref. na, where it was defined A C D — c' 6 H) 3 /v 2 (here we 

introduce the symbol c' 6 to distinguish it from c 6 in Eq. 0). with A 4 denoting the coefficient 
of the (WH) 2 operator in the SM Lagrangian. We fold c 3 = (1 + 5cg/2)/(l + 3c' 6 /2), 
which should be compared with the relation between c 3 and (valid at all orders in cq) 
appearing in Eq. (J 5 ]) . For small values of Cg one has A 4 ~ (m|/2u 2 ), hence c$ = c' 6 + 0(cg 2 ), 
so that the dependence of c 3 upon c' 6 is the same as in Eq. ([5]) up to small corrections. The 
parametrization of Eq. ([2]) is thus more convenient than that of Ref. [15] when Cq (or Cg) 
becomes large, since the formula for the trilincar coupling c 3 remains linear in cq. This is 
relevant for this analysis, since we anticipate that the results of Section [V] will constrain 
values of c 6 larger than 1 at the LHC. 


A. Modified power counting for the Higgs trilinear coupling 

The estimates of Eq. (|3]) shows that assuming the SILH power counting the modification 
of the Higgs trilinear coupling is expected to be of order (u//) 2 , i.e. of the same size as the 
shifts to other Higgs couplings. These latter however are already constrained from single- 
Higgs measurements to be close to their SM value, in particular the coupling of the Higgs to 
two vector bosons. It is thus interesting to ask whether there are scenarios, characterized by 
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a power counting different from the SILH one, where the Higgs trilincar coupling can get a 
large modification from NP effects while all the other Higgs couplings are close to their SM 
values, in agreement with the current LHC limits. Interestingly the answer to this question 
is affirmative: it is possible to imagine at least a few scenarios where the largest NP effects 
are in the trilincar coupling. 

A first possibility is a scenario where the Higgs is a generic bound state -i.e. not a Nambu- 
Goldstone boson- of some new strong dynamics. In this case no weak spurion suppression is 
required to generate 0 6 , and the naive estimate of c 6 is enhanced by a factor g 2 / A 4 compared 
to the SILH case, where we conveniently defined A 4 = m|/2u 2 ~ 0.13. The Higgs trilinear 
coupling (in SM units) is thus expected to be of order c 3 = 1 + 0(g 2 u 4 // 2 m 2 ), and large 
modifications are possible for g* large even if ( v/f ) 2 <C 1. For example, (v/ f) 2 = 0.05 gives 
ch ~ 0.05 and ~ 3.5 (g*/3) 2 , corresponding to c 3 ~ 0.93 + 3.5 (g*/3) 2 . The price to pay 
for this enhancement, however, is the tuning which is now required to keep the Higgs mass 
light, since naively one would expect m\ ~ m 2 . This is in addition to the 0(v 2 /f 2 ) tuning 
which must occur in the vacuum alignment even for a NG boson Higgs. Notice that similarly 
to Oq, other operators which previously required a breaking of the Goldstone symmetry to 
be generated will now be unsuppressed. In particular, c g and the coefficient of the operator 
WHB^B^ get enhanced by a factor (g 2 /X 2 ) compared to their SILH estimate, thus leading 
to modifications of the Higgs couplings to gluons and photons of order v 2 /f 2 times their 
(loop-induced) SM value. However, for (v/f) 2 -C 1 the current LHC constraints on these 
couplings can be easily satisfied. 

It is possible to avoid the tuning of the Higgs mass while keeping an enhanced trilin¬ 
ear coupling by envisaging a different scenario. Consider a theory where a new strongly- 
interacting sector, characterized by a mass scale m* and a coupling strength g*, couples to 
the SM sector only through the Higgs mass term portal: C mt = A HO. The Higgs field 
is thus elementary, and couples to the composite operator O, made of fields of the strong 
sector, with coupling A. This leads to the following estimates: 


AA4 ~ 



ch ~ 


A 2 u 2 

at P ’ 


A 2 v 2 y 2 _ A 3 v 2 

7.P 16^’ C6 ~9*V?’ 


( 6 ) 


where AA4 is the correction to the coefficient of the operator (H^H) 2 induced by the strong 
dynamics. The contribution to c u follows from a Higgs-dependent modification of the top 
quark kinetic term arising at the 1-loop level, and it is subleading compared to c#. Simi¬ 
lar effects are also generated for the lighter quarks proportional to their Yukawa coupling 
squared. Corrections to Eq. § of higher order in A are suppressed by powers of (A v 2 jrnf). 
An additional factor g 2 / 167 r 2 should be included in the above estimates if they arise from 
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diagrams involving a loop of strongly-coupled particles. 4 Since cq ~ (A/A4) ch, one can have 
an enhanced NP correction to the Higgs trilinear coupling for A > A4. However, requiring 
that no fine-tuning occurs in the coefficient of the quartic operator -i.e. in the physical 
Higgs mass- implies A < \f\ 4 <7*. Furthermore, only for (Au 2 /m*) < 1 one can make sense 
of the series in powers of the Higgs held, since when this inequality is saturated all orders in 
A become equally important and perturbation theory is lost. From the above considerations 
it follows the bound c 6 < 0(1). In this scenario it is thus possible to avoid tuning the Higgs 
mass but the corrections to the Higgs trilinear are at most of 0(1); larger enhancements 
require fine-tuning. For example, <7* = 47 r, m* = 500 GeV and A = ~ 4.5 gives 

ch ~ 0.03 and c 6 rvj t 

These examples show that it is possible, in scenarios characterized by a power counting 
different from the SILH one, to have large modifications to the trilinear coupling while 
keeping the other Higgs couplings close to their SM values. In some cases, like that of a 
generic composite Higgs, c$ can be of order 1 or larger without invalidating the perturbative 
expansion in powers of the Higgs held. 

B. The relevance of higher-order operators 

Let us now investigate more quantitatively the validity of our effective Lagrangians and 
discuss the importance of higher-order operators. We will do so by adopting the SILH power 
counting. We already mentioned that in the case of the linear Lagrangian, for a given process, 
operators with higher powers of the Higgs doublet imply corrections suppressed by additional 
factors (■ v/f ). This means that if v/f ~ 1 the linear Lagrangian cannot be consistently 
used anymore, although we can still rely on the non-linear one. As regards the derivative 
expansion, the same considerations apply instead to both descriptions. The perturbative 
parameter in that case is ( E/m *), so that naively higher-order operators become important 
at E ~ m*, where the effective description breaks down. However, in certain processes it 
can happen that the leading operators are suppressed due to symmetry reasons, and the 
higher-order ones become relevant at large energies below the cutoff scale. This is in fact 
the case of double Higgs production via gluon fusion, where the dimension-6 operator O g is 
suppressed by two powers of a weak spurion if the Higgs is a NG boson. 

4 This is what happens, for example, when the strong dynamics consists of a single real scalar field S, 
neutral under the SM gauge group. By requiring invariance under the parity S —► —S, the Lagrangian 
describing the strong dynamics is C = (d^S) 2 /2 — m 2 s S 2 /2 — AsS 4 /4, and O = S 2 . One can thus identify 
m* = ms and g* = vWs- 
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The dimension-8 operators relevant for double Higgs production via gluon fusion are: 


A£ 8 d 


m 


w L 


c gD o (D p H^DPH)G a ^G 


a fiu 


+ c gD 2 {rTDprfDPH - AD^H^D u H)G; a G a u a 


( 7 ) 


They mediate scatterings of two gluons into two Higgs bosons with total angular momentum 
equal to respectively 0 and 2, and can therefore lead to different kinematic configurations, 
as we will discuss later on. Both operators, however, are similar from the power-counting 
viewpoint. Neither of the two breaks the shift symmetry of a NGB Higgs, and thus no 
spurion suppression appears in the estimate of their coefficients: 5 


CgD 0,2 





( 8 ) 


We can now compare the contributions of the various operators to double Higgs produc¬ 
tion via gluon fusion. 6 7 By using the estimates of Eqs. the scattering amplitude can 

be schematically expressed as 


A (gg hh) ~ (^) x 


P 


y 2 t[^ + 0[-\\+gl{E)+gl{E) + 


where 


2/r ^ _ 4 tt E 2 A 2 E 2 

gpE) ~ Cg - - ~ —— 

«2 v z mi 


g%{E) ~ CgD 0,2 — 


4tt £ 4 


«2 


glE 4 


mz 


(9) 


( 10 ) 


The first term in the square parentheses of Eq. (|9]) corresponds to the SM contribution 
plus the 0(y 2 / f 2 ) correction implied by Oh and O u . 1 The second and third terms in the 
parentheses denote the contributions of respectively O g and 0 9 do , 0 9 d 2 , and thus define the 
coupling strengths of the interactions mediated by these operators. Since these couplings 
grow with E, at sufficiently large energies, yet below the cutoff scale, they can become larger 
than the top Yukawa coupling. It is thus possible to obtain NP corrections larger than the 
SM contribution within the validity of the effective Lagrangian. This should be contrasted 
with the corrections from dimension-6 operators to the on-shell production and decay rates 
of the Higgs boson: those are at most of order 0(u 2 // 2 ), hence always smaller than the SM 


5 This estimate assumes that O g n 0 and O g n 2 are generated at 1-loop by the New Physics. The same 
assumption has been done on O g in the estimate of Eq. (|3j). 

6 R.C. would like to thank Riccardo Rattazzi for illuminating discussions on the validity of the effective 
theory in scattering processes. For a related analysis of WW scattering, see: R. Rattazzi, talk at BSM 
physics opportunities at lOOTeV , CERN, February 2014. 

7 Here we neglect the logarithmic energy growth of the triangle diagram with the tthh quartic interaction 
generated by O u . The correction implied by Oq enters through the SM triangle diagram and is further 
suppressed at high energies, see next section. 












term. We thus see the potential advantage of studying 2 —)• 2 scatterings at high energies 
compared to the single-Higgs measurements performed during Runl at the LHC: going off- 
shell at higher energies one can directly probe the strength of the New Physics dynamics 
in the regime in which it gives large effects [TH]. In practice, as we will discuss more in 
detail in the following, such regime is difficult to exploit in the process gg —>■ hh —>■ 6677 
considered in this work. This is due to the steep fall off of the gluon pdfs at large x and 
to the small final-state branching ratio, which strongly limit the exploration of the region 
with large hh invariant mass. For A > ^Jg*yt the coupling ge(E) is the first to become larger 
than y t at a scale E ~ y t f (g*/ A). The smaller A is, the higher the crossover scale becomes, 
until E ~ y t f ( g*/yt ) 1//2 is reached where g%{E) ~ y t . Above this energy, for A < y/g^Ut, the 
largest effects come from the dimension -8 operators. The validity of the effective Lagrangian 
extends up to energies E ~ m*, so that the coupling strengths are bounded by 

9e(E)<\<g*, g 8 (E)<g *. ( 11 ) 


As we anticipated, studying the kinematic region where higher-order operators give a 
contribution larger than the SM one is challenging in double Higgs production via gluon 
fusion. However, it is still interesting, and relevant for the analysis carried in this work, to 
ask at which scale the dimension -8 operators become more important than the dimension -6 
ones, independently of their absolute size. From Eq. (10) it follows that g§(E) 9%{E) for 

E ~ A/; at this scale g$(E) ~ A 2 /g*. Hence, a further condition to be satisfied in an analysis 
which includes only dimension -6 operators neglecting those with dimension 8 is: 


Se (E) < - . 
g* 


( 12 ) 


Let us then indicate with 6 the precision obtained in such analysis on C 2 g = c g (An / 0 . 2 ), by 
making use of events with invariant masses up to E. This means that the smallest value 
of ga probed is g 

171171 VS(E/v), so that Eqs. (|ll| , (Il2j) can be re-cast into the following 

constraints: 


A grain (13) 

A 2 

grain • (14) 

g* 

These inequalities define the region in the (g*, A) parameter space, sketched in Fig. [TJ which 
can be sensibly probed, within the validity of the effective theory, by an analysis including 
only dimension -6 operators. Clearly, smaller values of g m i n , obtained by increasing the 
precision 5 for a fixed energy E, lead to a larger viable region. Notice that the two cases 
A — yt (fully composite £r) and A = g*y t (partially composite ij, and £r) can be probed 
only if g min < y t . This is not the case of course if the analysis does not have enough 
sensitivity to probe the SM cross section. 
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FIG. 1: Cartoon of the region in the plane ( g *, A /g*), defined by Eqs. (13),(14), that can be probed 
by an analysis including only dimension-6 operators (in white). No sensible effective field theory 
description is possible in the gray area (A < g m in)i while exploration of the light blue region 
{dmin < A < yjg*gmin) requires including the dimension-8 operators. 



FIG. 2: Feyman diagrams contributing to double Higgs production via gluon fusion (an additional 
contribution comes from the crossing of the box diagram). The last diagram on the first line 
contains the tthh coupling, while those in the second line involve contact interactions between the 
Higgs and the gluons denoted with a cross. 


C. Cross section of double Higgs production 

We can now discuss our parametrization of the cross section of double Higgs production 
via gluon fusion. We will use the non-linear Lagrangian Q and start by neglecting higher- 
derivative terms (which correspond to dimension-8 operators in the limit of linearly-realized 
EW symmetry). The effect of the neglected derivative operators will be then studied by 
analyzing their impact on angular differential distributions and shown to be small in our 
case due to the limited sensitivity on the high m^h region. 

The Feynman diagrams that contribute to the gg —» hh process are shown in Fig. [2] Each 
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diagram is characterized by a different scaling at large energies \/I = rrihh 3> mt, We 


find: 



/ mf 
I log -y + *7T 




2 


(15) 




where Ma, Ma are the amplitudes of respectively the box and triangle diagram with Higgs 
exchange, AauI is the amplitude of the triangle diagram with the tthh interaction, and 
As and A 4 are the amplitudes of the diagrams with the Higgs-gluon contact interactions. 
One can see that each NP contribution affects the rrihh distribution in a different way. In 
particular, the diagrams that depend on the Higgs trilinear coupling c 3 are always suppressed 
at large s, and their contribution affects the process mostly at threshold. Modified values 
of the top Yukawa coupling c t and the non-linear interactions C 2 t and C 2 g , instead, tend to 
increase the cross section at higher invariant masses. Finally, including the dimension-8 
operators would lead to an additional contribution to A 4 growing as s 2 and distort the tail 
of the rrihh distribution. A shape analysis can thus help to differentiate the different effects 
and break the degeneracy of the total cross section on the Higgs couplings. This will be our 
strategy in the study of double Higgs production discussed in the next section, where we 
will use mhh as the main kinematic variable to characterize signal events. 

By focussing on gluon fusion, the total cross section for the process pp —> hh can be 
written as a simple polynomial of the parameters of the effective Lagrangian: 8 



The LO value of the numerical coefficients Ai and of the SM cross section &sm is reported 


8 The contribution from Vector Boson Fusion is smaller by at most one order of magnitude, and can be 
isolated by selecting the number of jets in the final state. See for example Refs. mm for up-to-date 
studies with mu = 125 GeV. 
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VG 

o+m [fb] 

Ai 

a 2 

A 3 

a 4 

A 5 

Aq 

a 7 

14 TeV 

16.2 

2.13 

10.1 

0.300 

21.8 

188 

-8.62 

-1.43 

100 TeV 

874 

1.95 

11.2 

0.229 

16.0 

386 

-8.32 

-1.18 

V* 

As 

Ag 

A 10 

An 

A 12 

A 13 

All 

A 15 

14 TeV 

2.93 

21.0 

59.8 

-9.93 

-23.1 

4.87 

10.5 

96.6 

100 TeV 

2.55 

16.9 

52.4 

-7.49 

-17.3 

3.55 

8.46 

87.5 


TABLE I: Coefficients of the fit of the total pp —> hh cross section via gluon fusion given in Eq. (16). 
They have been computed at LO for the 14 TeV LHC and for a future collider with yfs = 100 TeV. 


in Table [T] for hadron colliders with center-of-mass energy a/s = 14 TeV and 100 TeV. 9 
They were computed with our dedicated C++ code linked to QCDLoop [21] and to the 
LHAPDF routines [22], by setting rr+ = 125 GeV, rn t = 173 GeV and using the CTEQ611 
parton distribution functions. The factorization and renormalization scales have been fixed 
to Q = VI = m hh . 

Besides mhh, the other kinematic variable that characterizes the gg -+ hh events is the 
angle 6 between either of the two Higgses and the beam axis in the center-of-mass frame. 
By total angular momentum conservation, the scattering amplitude can be decomposed into 
two terms, M 0 and M 2 , describing transitions with respectively J z = 0 and J z = ±2, where 
J z is the projection of J along the beam axis. The amplitude M 2 receives a contribution only 
from the box diagram and from the dimension-8 operator O g m (through the last diagram 
of Fig. [2]). All the other diagrams, instead, mediate J z = 0 transitions. The explicit 
expression of M 0 and M 2 is reported in Appendix [A] for convenience. While transitions with 
all integer values of the total angular momentum J can occur, the leading contributions to 
M 0 and M 2 come respectively from J = 0 (s-wave) and J = 2 (d-wave). A simple angular 
momentum decomposition then shows that, up to small corrections, the angular dependence 
is M 0 ~ const., M 2 ~ sin 2 9. 10 In the SM, the contribution of M 2 to the cross section is 
always small, as illustrated by Fig. [3] it is negligible at the peak of the mhh distribution 
(' m hh 400 GeV) and smaller than ~ 20% on its tail (mhh ~ 700 GeV). A shift in the top 
Yukawa coupling due to New Physics modifies the value of the box diagram, but cannot 


9 We computed the Aj’s by performing a corresponding number of MC simulations, each with ~ 10 6 events. 
We estimate the statistical uncertainty on the A,;’s to be at the 10 “ 2 level. This is less accurate than one 
could naively deduce due to cancellations that occur when extracting some of the coefficients, but still 
much smaller than the theoretical error from PDFs and higher-order QCD corrections (see footnote 20 and 
Ref. 0 ). When we consider different mhh categories later in the analysis, the fit is individually performed 
in each mhh bin through a similar procedure, thus maintaining the same level of statistical uncertainty. 

10 The angular dependence of a scattering amplitude with J = j and J z = m is given by the Wigner function 


12 







1 

S' 

CO 

o 

o 

;o 

o -1 

S 10 


10 2 


0 0.2 0.4 0.6 0.8 1 

cose 

FIG. 3: Differential cross sections obtained by including only the contribution of Mq (dashed blue 
curve) or M 2 (dotted orange curve) in the SM, as functions of cos 6. Both curves are normalized 
to the total SM cross section. The partonic center-of-mass energy has been fixed to yfi = 400 GeV 
in the left plot and to y/I = 700 GeV in the right plot. 
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change the above conclusion (unless extreme shifts are considered). It is thus interesting to 
ask whether the dimension -8 operator O g £> 2 can give a sizable contribution to M 2 through 
the last diagram of Fig. [2j If this were the case, one could reach a better sensitivity on O g o 2 
by performing a suitable analysis of the angular distributions of the Higgs decay products. 
Unfortunately, we find that the effect is numerically small. This is illustrated in Fig. [4| 
where we show the isocurves of r = cr g r, 2 /a tot in the plane [rrihh: cos 9 min ). We defined r to 
be the ratio of the cross section induced by O g n 2 alone, a g o 2 , to the total cross section, cr tot) 
obtained by adding the contribution of O g o 2 to the SM. Both cross sections are computed 
at the partonic level for the process gg —» hh and by integrating over angles 9 min < 9 < 7 t/2 . 
We set the coefficient c g n 2 equal to its NDA estimate ([ 8 ]) and choose as benchmark values 
/ = 635 GeV, m* = 1.9 TeV, which correspond to g* = 3, £ = ( v/f ) 2 = 0.15. For r ~ 0.5 
the contribution of O g D 2 is as large as the SM one. The plot of Fig. [4] shows that this occurs 
at the crossover scale rrihh ~ 1.3 TeV, in agreement with the naive estimate of section IIB 
rrihh ~ Utf (g*/yt) 1/2 - Increasing the value of 9 mm (i.e. decreasing cos 9 mm ) enhances the 
importance of the J z = 2 component of the cross section, hence that of O g D 2 , but the effect 
is never large. For example, the value of the crossover scale is reduced at most by ~ 10%. 
We thus conclude that, although an analysis of angular distributions could in principle help 
disentangling the effects of dimension -8 operators, in practice there is little leverage, and the 
efficacy of such strategy is further reduced in the case of the 6677 final state by the limited 
range in rrihh which can be realistically probed. In our analysis of gg -7- hh —y 6677 we will 
thus make use of the rrihh distribution as the main tool to probe the effects induced by New 
Physics, neglecting for simplicity angular distributions. 
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FIG. 4: Isocontours of the ratio r = (JgDijvtot , defined in the text, in the plane (mhh, cos 9 m i n ). 


III. DOUBLE HIGGS PHENOMENOLOGY AT 14TEV AND 100 TEV 

The presence of various Higgs decay channels with a non-negligible branching ratio allows 
the exploration of the Higgs properties in several final states. This is especially true in single 
Higgs processes, where the relatively large production cross section compensates for the small 
decay probability in some of the cleanest channels, such as h —>■ 77 and h — > ZZ* —>■ 4 l. In 
the case of double Higgs production, however, the low signal rate forces one to consider only 
final states with a sizable branching ratio. This is particularly so at the 14 TeV LHC, where 
the NLO total SM production rate is around 37 fb -1 , but remains true even at a hypothetical 
future 100 TeV machine, where the enhanced rates are unavoidably accompanied by larger 
backgrounds. In order to obtain a large enough branching ratio, at least one of the two 
Higgses must decay into a bb pair. For the second Higgs different choices seem possible 
and have been considered in the literature, namely h —s- bb, h —> WW*, h — > r + r~ and 
h — > 77 . The channel hh —>■ bbbb has the highest signal rate ( BR ~ 33.3% in the SM), 
but is plagued by a large QCD background. Even imposing four 6-tags, it seems hard to 
exploit and could require a sophisticated search strategy mum- The channel with the 
second highest rate is hh —> bbWW*, with a branching ratio BR ~ 24.9% in the SM. Its 
observation is also threatened by a large background, mainly coming from ti EE]. It has 
been recently claimed that a good signal-to-background ratio can be obtained by using jet 
substructure techniques and focusing on a very specific region of the parameter space where 
the two Higgses and their decay products are highly boosted [6]. Although encouraging, the 
results of this analysis suggest that the bbWW* final state can be observed at the 14 TeV 
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LHC only with its high-luminosity extension. The other channel which has been extensively 
studied in the literature is hh —> bbr + T~ ’ PI5U7U9UI3]. It is very promising and potentially 
relevant for the 14TeV LHC, since its SM branching ratio is sizable, BR ~ 7.35%, and good 
signal-to-background ratios seem to be achievable while keeping a relatively large number 
of signal events. Its actual significance relies however on the ability to reconstruct the tau 
pair, and will have to be fully assessed by an experimental study. 

For the purpose of our study we will focus on the channel hh —> 77 bb (BR ~ 0.264% in 
the SM), which has been considered to be the cleanest one despite its small rate. As shown 
by previous theoretical studies 00 IEI flO] as well as a recent non-resonant experimental 
search at /~s = 8 TeV [23] and a study for the /s = 14TeV high-luminosity LHC [24] 
by ATLAS, 11 the analysis strategy to observe this decay mode is relatively simple and 
straightforward. This allows us to avoid unnecessary complications and to concentrate 
primarily on the interpretation of the search in the context of the effective held theory 
description. Clearly, exploring the other available channels from a similar perspective could 
improve significantly the sensitivity. We leave this interesting follow-up for a future study. 

Before proceeding with the details of our analysis it is useful to briefly summarize the 
properties of the kinematic distributions and discuss the main differences between the 14 TeV 
LHC and a future 100 TeV collider. 

It is well known (see for example Ref. [12] ) that in the SM a cancellation between the box 
diagram and the triangle diagram involving the Higgs trilinear coupling leads to a depletion 
of the signal at threshold (rrihh > 250 GeV). The peak of the distribution is essentially 
determined by the fast decrease of the gluon parton distribution functions and is located 
around rrihh — 400 GeV. This conclusion is valid independently of the collider center-of-mass 
energy /~s . The main difference between the LHC and a higher-energy collider is a rescaling 
of the overall cross section (the SM cross section at a 100 TeV machine is around 40 times 
bigger than the one at the 14TeV LHC), with small effects on the rrihh distribution. As it 
can be seen from Fig. [5j this latter is modified significantly only on its tail (rrihh > 700 GeV), 
which is enhanced at larger collider energies due to the higher luminosity of gluons. In fact, 


as we will discuss in section IV A an important change is also present in the pseudo-rapidity 
distribution of the Higgs bosons, which is related to the amount of boost determined by the 
initial parton energies. 

Due to the above considerations, one would naively expect that an analysis similar to the 
one designed for 14 TeV would continue to work well at 100 TeV. This is in fact only partially 
true. Indeed, the enhanced signal rates can impact the analysis results in two ways. First, 


11 See Ref. [25] and Refs. [231 f2S] for respectively theoretical and experimental analyses of the resonant 
di-Higgs production in the 6677 final state. 
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FIG. 5: Left plot: Normalized differential cross section for pp —> hh in the SM as a function of 
the invariant mass of the two Higgs bosons. The solid and dotted lines correspond respectively to 
y/s = 14 and 100 TeV. Right plot: Same as on the left but with logarithmic scale. 


they will improve the sensitivity on the parameters, such as the Higgs trilincar coupling, 
by simply reducing the statistical uncertainty even if the analysis strategy is not modified. 
Second, due to the larger cross section and enhanced tail, much higher values in rrihh will be 
accessible, which are potentially more sensitive to higher-order operators growing with the 
energy. For this reason, fully exploiting the potential of a 100 TeV collider requires some 
modification of our 14 TeV analysis strategy in order to better reconstruct events with a 
higher boost. Although we will not present a fully optimized analysis for the 100 TeV, we 
will give a first assessment of how much our final results can improve when jet substructure 
techniques are used. 

As a first step in this direction, it is useful to get an idea of the reach that can be 
achieved on rrihh and pr(h) at different colliders and for different search channels. A complete 
assessment of this point would require specifying completely the scenario we are interested 
in and the size of the possible new physics contributions to the cross section. To get a rough 
approximation, however, we can estimate the reach by demanding that a few events (we 
choose 5 for the estimate) are still present in the tails of the SM distributions above the 
considered rrihh (or PtO 1 )) value. To be more realistic we also assume a 10% overall signal 
efficiency due to kinematic cuts. The estimates can be easily extracted from the plots of 
Fig. [6j which show the integrated differential cross sections for the SM signal pp —>■ hh as 
a function of the lower cut on rrihh and of the minimum pr{h). By including the branching 
ratio of the most important decay modes we obtain the results shown in Table [TT| 
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FIG. 6: Integrated differential cross sections for pp —> hh in the SM as a function of the minimum 
hh invariant mass, (m,hh)m in (left plot), and the minimum Higgs transverse momentum, (pr(h)) m \n 
(right plot). The upper (lower) curve is for -y/s = 100TeV (14 TeV). The cross sections have been 
computed at NNLO using the k-factors of Eq. (EDI¬ 


channel 

bbbb (33.3%) 

bbWW* (24.9%) 

bbr + r (7.35%) 

7766 (0.264%) 

Cross section 
m hh [GeV] 
p T (h) [GeV] 

> 0.05 fb 

< 1340 (4290) 

< 575 (2000) 

> 0.067 fb 

< 1280 (4170) 

< 550 (1890) 

> 0.227 fb 

< 1039 (3235) 

< 440 (1430) 

> 6.31 fir 

< 558 (1552) 

< 210 (664) 


TABLE II: Lower limits on the pp —>• hh cross section in the SM and upper limits of the phase 
space that guarantee 5 signal events. We assume an integrated luminosity L = 3000 fb -1 at the 
14 TeV LHC and a 10% efficiency due to kinematic cuts. The numbers in parenthesis refer to a 
100 TeV collider, with the same integrated luminosity. The cross sections have been computed at 
NNLO using the k-factors of Eq. ( [17] ). 

The reach on is rather limited in the 6677 channel due to the small signal rate, 
and even at a 100 TeV collider it does not extend much beyond ~ 1.5 TeV. Other channels 
with larger rates will be crucial to push further the exploration of higher invariant masses. 
Similar considerations apply also for the reach on pr(h)- The maximal value of pr(h) can 
be used to estimate whether a given search can benefit from jet substructure techniques or 
a simple jet analysis is enough. The angular separation between two partons coming from 
the Higgs decay scales like A R ~ 2 mh/prih)- In order to resolve the jets with the standard 
techniques we must demand A R > Ratlas,cms, where Ratlas = 0.4 and Rcms = 0.5 
are the reconstruction cones used by the ATLAS and CMS collaborations respectively. It 
is easy to see that for p T > 500 — 600 GeV the two partons are not often resolved as two 
separated objects and jet substructure techniques are likely to be helpful. The results in 
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Table [TT] show that the size of the boosted phase space at 14 TeV is almost negligible for the 
77 bb channel and very limited for bbbb , yylVIV* and bbr + r ~. The story changes however 
for an 100 TeV collider. In this case the SM signal can be probed up to rrihh ~ 1.5 TeV in 
the 77 bb channel and up to ~ 4 TeV for bbbb. In these kinematic regions many boosted 
Higgses are produced and jet substructure techniques are crucial to reconstruct them. We 


will analyze this aspect in more details in Section IV A 


IV. ANALYSIS OF pp ->• hh 77 bb 


After discussing the general properties of the pp —>■ gg —> hh process, we now focus on the 
6677 final state and describe our analysis strategy. We generated the signal events through 
our own dedicated code described in section II C| The code includes the 1-loop diagrams of 
Fig. | and is based on the parametrization of Eq. Q. It is available from the authors upon 
request. By using the CTEQ611 PDF’s (LO PDF with LO a s ), setting the factorization and 
renormalization scales to Q = mhh and rrih = 125 GeV, mt = 173 GeV, we find the following 
LO cross sections for the Standard Model: 16.2 fb and 873.6 fb for respectively 14 TeV and 
100 TeV collider energies. In order to (partially) include the NLO and NNLO corrections, 
we rescale the cross sections by the k-faetors 


^14TeV — 2.27 , &100TeV — 1-75 (17) 

which were computed for the SM in Ref. EH . 12 

The main backgrounds that we considered are the non-resonant processes 6677 and 77 jj 
and the resonant processes bbh, Zh and tth (with the subsequent decays h —> 77 , Z —> bb, 
and tt —>■ bb -\- X). Further backgrounds involving fake photons from jets have been ne¬ 
glected as their estimate is beyond the scope of this work. Experimental analyses of single 
Higgs production have shown that similar processes in that case can be safely reduced to a 
subleading level, and we thus assume that this will be possible for double Higgs production 
as well. We generated all the backgrounds with MADGRAPH5_aMC@NLO v2.1.1 [33J by 


12 Compared to the NLO calculation (performed in Ref. [21]), the NNLO corrections give a ~ 20% enhance¬ 
ment of the total cross section. Notice that both the NLO and NNLO corrections to the total rate have 
been computed in the infinite top mass approximation. It is well known that at LO this approximation 
leads to very inaccurate kinematic distributions, in particular the mhh one. For this reason the k-factors 
used in our analysis must be considered only as a crude approximation. Recently, a step forward towards 
a fully differential NLO calculation has been done in Ref. [2H], where real emission diagrams were com¬ 
puted at 1-loop, while the finite part of the two-loop virtual corrections was extracted by resorting to the 
infinite top mass approximation. For an estimate of the finite-top mass effects and of the dependence of 
the k-factor on a cut on -/§, see Refs. mm- See also Ref. [35] for the resummation of threshold effects 
in the SM. 
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switching off virtual corrections (i.e. working in LO mode). The output has been inter¬ 
faced with PYTHIA v6.4 [M] for parton showering and hadronization, and with FastJet 
v3.0.6 [3Sj for jet clustering. The factorization and renormalization scales have been set to 
the default dynamic scale in MadGraph5. Further details about the generation can be 


found in Appendix B 1 The backgrounds 6677 , bbh and Zh have been matched up to one 
additional jet via the kx -jet MLM matching [[36] to partially account for NLO effects. In 
the case of 6677 we found that allowing for the extra jet increases the cross section by a 
factor ~ 2. We found that this is in good agreement with a complete NLO computation per¬ 
formed with MADGRAPH5_aMC@NLO, indicating that real emissions in this case represent 
the bulk of the NLO correction. Considering that 6677 is the dominant background after all 
cuts, including this effect is very important, but it has been omitted in most of the previous 
studies. More details about our generation of 6677 and a thorough discussion of the NLO 


correction are given in Appendix |B 2[ For the resonant tth background we apply instead 
a simple k-factor rescaling of the LO cross section to the NLO value reported in Ref. m- 
Finally, the 77 jj background was simulated without matching and a k-factor fc 77? -j = 2 was 
applied (this is a more conservative choice than the one proposed in [33] based on an NLO 
simulation). 

To estimate the sensitivity of our analysis we should, ideally, take into account parton- 
shower/hadronization and smearing in every point of our parameter space. This procedure, 
however, would be computationally too expensive. Instead we adopt the following simplified 
approach. We fully extract the signal rate after cuts only for the SM point by performing a 
hadron-level analysis. For the other points in the parameter space, we apply the same type 
of analysis to the parton-level samples and we rescale the signal by the hadron-to-parton 
cross section ratio computed for the SM: 

r had \ 

(18) 


had I ^ part I 

BSM| w / cuts BSM| w / cuts 


cr, 


SM \ 
part J 
°SM / 


w/ cuts 


We use the fact that such ratio is approximately the same for both the SM and in a generic 
point of the BSM parameter space. The rescaling is performed individually for each of the 
bins in rrihh of Tables [Y] |VII| and [VIII[ Our analysis takes correctly into account possible non- 
universal signal effects coming from distortions of the kinematic distributions at the partonic 
level. The rescaling then approximately includes the effect of showering and hadronization. 
For the jet substructure analysis of Section IV A , on the other hand, we follow a different 
procedure, since there is no parton-level cross section we can use after cuts. We thus compute 
the BSM signal by starting with the SM one computed at the hadron level and multiply by 
the ratio of BSM over SM cross sections before cuts at the parton level 


had I ^ had I 

BSM| w / cu t s SM | w / cu t s 


X 


a 


a, 


part 

BSM 

part 

SM 


w/o cuts 


(boosted analysis). 


(19) 
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This ratio is expected to be approximately the same at the partonic level before cuts and at 
the hadronic level after cuts if one considers sufficiently narrow bins in rrihh, since to very 
good accuracy the latter is the only variable that controls the kinematics of the signal. We 
checked that this procedure is accurate at ~ a few % level in the high-m^ categories used 
for the boosted analysis. 

We designed different strategies for 14 TeV and 100 TeV colliders. We begin by describing 
the one at 14 TeV. Events are triggered by demanding exactly two isolated photons satisfying 
the minimal reconstruction requirements 


Pri'y) > 25 GeV, |p( 7 )|<2.5. (20) 

To forbid extraneous leptonic activity we veto events with isolated leptons (electrons or 
muons) satisfying 

Pt(1) > 20 GeV, \r)(l)\ <2.5. (21) 

Similar isolation criteria are applied to both photons and leptons. A photon (lepton) is 
considered isolated if the surrounding hadronic activity within a cone of size R = 0.4 (0.3) 
satisfies Pt( 7 )/(pt( 7 ) +pr(cone)) > 0.9 (0.85). We have checked that more sophisticated 
photon isolation criteria, like the one proposed in [38], give similar results. We then cluster 
the events into R = 0.5 anti -k T jets (39]. We accept only jets with 

Pt(J)> 25 GeV, \ V (j)\ < 2.5 (22) 


and require that at least two of them are 6-tagged (the leading two 6-jets are selected if 
more than two 6-jets exist). We assume 70% efficiency for 6 tagging e& = 0.7, corresponding 
to 1% of mis-tag rate = 0.01 14011421 . 13 and 80% efficiency for photon tagging e 7 = 
0.8 [121 S3] 14 . After this step, an event consists of two isolated photons, two 6-tagged jets, 
and possible additional jets (whether or not 6-tagged). 

Once the reconstruction of the basic objects by the above procedure is done, we apply 
the set of cuts described in the following. We first restrict the events to those with two hard 
photons and two 6-tagged jets satisfying 


Pt>(&),Pt>( 7 ) > 50 GeV, Pt<(&),Pt<( 7 ) > 30 GeV, 
60 < m2 co < 200 GeV, 60 < m™° < 200 GeV, 


(23) 


13 More specifically, we tag with 70% probability only those jets which contain a 6-liadron with transverse 
momentum larger than 5G eV/Rj e t- For Rj e t = 0.5 this translates into a minimum 10 GeV transverse 
momentum. When the jet substructure is used (see section IV A), Rj et is set to A R(bb), corresponding 
to the distance between two leading subjets after BDRS declustering. 

14 Typical rejection rates for fake photons from jets at this working point are in the range 0.1 — 0.5%. We 
will not use this number however, since backgrounds from fake photons are not included in this analysis. 
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Signal (SM), 14 TeV 77 bb, 14TeV tth, 14 TeV 



0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 

A R(b,b) A R(b,b) A R(b,b) 


FIG. 7: Angular distributions (in arbitary units) after the cuts of Eq. (23) of the SM signal and 
the two dominant backgrounds, 77 bb and tth , for s/s = 14 TeV. The dashed lines indicate to the 
cuts of Eq. (f24l) . 


where Pt>{i) (Pt<{ 7 )) denote the transverse momentum of the hardest (softest) photon, 
and pr>(b), pr<(b) are similarly defined for the two 6-jets. At this stage, the broad mass 


windows in Eq. (23) are placed merely to allow for a fair comparison of the signal and 


background cross sections. 

In the signal, the bb and 77 subsystems tend to be in a back-to-back configuration, 
and the angular distance between two photons (and similarly between two 6 -jets) is of order 
~ 2 mh/pr{h) ~ C?(l) in the majority of the phase space. While the 77 subsystem in resonant 
backgrounds has a kinematics similar to the signal, the different origin of the bb pair in each 
process can be used to distinguish them from the signal. The angular distributions of the 
signal and of the two dominant backgrounds are shown in Fig. [7j We find that the following 
cuts (indicated by the dashed lines in Fig. [7j) 


A R(b,b) < 2 , Ai 2 ( 7 , 7 ) < 2 , AR(b, 7 ) > 1.5, 


(24) 


can efficiently reduce the background, especially 6677 , while retaining most of the signal. 15 
As a final cut, we restrict the invariant masses of two 6 -jets and of the two photons to the 


15 At leading order the 6677 process is initiated by qq and gg and proceeds through diagrams where the two 
6 ’s tend to emerge with a large relative angle. Diagrams initiated by gq and gq at next-to-leading order 
lead to topologies which can more easily fake the angular configuration of the signal but account only for 
a minor fraction of events. 
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Vs = 14 TeV 

hh 

6677 77 jj tth bbh Zh 

After selection cuts of Eq 

After A R cuts of Eq. (24 

After m[f co cuts of Eq. ( 

00 ,77 ^ 

• ( 

) 

25 

23 

) 

) 

25.8 

17.8 

12.8 

6919 684 130 7.2 25.4 

1274 104 29 1.2 15.8 

24.2 2.21 9.9 0.40 0.41 


TABLE III: Cut flow for the SM signal and the various backgrounds at \/s = 14TeV. The values 
correspond to the number of events, assuming an integrated luminosity L = 3abW All recon¬ 
struction efficiencies and branching ratios to the final state under study are included. 



rsM [ab] 

A 1 

d -2 

A 3 

A4 

A 5 

Afi 

h 7 

14 TeV 

4.28 

1.70 

10.7 

0.117 

6.11 

217 

-7.56 

-0.819 

100 TeV 

92.9 

1.59 

12.8 

0.090 

5.20 

358 

-7.66 

-0.681 

VS 

d -8 

A9 

A10 

An 

A 12 

A 13 

An 

A15 

14 TeV 

1.95 

10.9 

51.6 

-3.86 

-12.5 

1.46 

5.49 

58.4 

100 TeV 

1.83 

9.25 

51.2 

-2.61 

-7.35 

1.03 

4.65 

65.5 


TABLE IV: Coefficients of the fit of the total signal rate (r = a x BR(hh —> 6677 )) obtained 
after all cuts at 14TeV and 100 TeV. The fit is based on a parametrization analog to that given in 
Eq. (16) for the cross section. The SM rates rsM include the NLO k-factors of Eq. ©■ 


Higgs mass window 16 


105 < m^ co < 145 GeV , 120 < m™° < 130 GeV. 


(25) 


The resulting cut flow is shown in Table |III| We report in Table |IV| a fit of the signal rate 
(r = a x BR{hh —)■ 6677 )) obtained after all cuts based on a parametrization analog to that 


given in Eq. (16) for the cross section. 


With our analysis strategy, 77 bb and tth are the two major backgrounds. The latter tends 
to produce extra jets from the hadronic decay of the W’s. One could thus consider applying 
a veto on the extra hadronic activity to enhance the signal significance. The potential impact 
of a jet veto can be seen in Fig. | 8 j For example, further restricting the events to the region 


16 The width of the interval 120 < m!V 0 < 130 GeV corresponds to ~ 3 times the experimental resolution 
on photon pairs, see [MlHS]. The same mass window was adopted by the recent CMS study of di-Higgs 
resonant production Jfjjj], see also Ref. 33|. The width of the interval 105 < m 1 '| co < 145 GeV corresponds 
to ~ 2 times the experimental resolution on 6 -jet pairs from Higgs decays (after correcting for resolution 
effects), see [Ml IS]- We do not smear photons in our simulation, nor include a specific efficiency for the 
reconstruction of the photon pair. The mass window on m^° is on the other hand sufficiently wide that 
almost all of the signal is retained, so that the efficiency would be close to 100 % even including a finite 


energy resolution on photons. The efficiency of the cuts of Eq. (25) reported in Table III corresponds, in 


the case of the signal, to the efficiency for the reconstruction of the bb pair (equal to ~ 0.72). 
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14 TeV | | 14 TeV | 




FIG. 8: Multiplicity of jets (left plot) and of reconstructed hadronic W bosons (right plot) after all 
cuts at \fs = 14 TeV. The solid, dotted and dot-dashed curves denote respectively the SM signal, 
77 bb and tth. 


with N(jets) < 4, in addition to all the previous cuts, can remove roughly 80% of the tth 
background while keeping ~ 70% of the signal. Alternatively, one could look for hadronic 
IT’s by iteratively forming jet pairs with invariant mass lying in a given window around the 
W mass, for example we will use (70,100) GeV in the following. As illustrated in Fig. [8j 
vetoing hadronic W 's in addition to our cuts can reduce ~ 50% of the tth background while 
retaining more than 90% of the signal. We find that applying either of these cuts leads to 
a modest increase in the significance, at the cost of reducing the number of signal events. 
We thus decided not to exploit any form of extra-jet vetoing at 14 TeV, motivated by the 
necessity of retaining as many signal events as possible. 


As discussed in Section IIC , several diagrams with different energy scalings contribute to 
the signal. We can thus use the invariant mass of the reconstructed hh system to differentiate 
the various effects and improve the sensitivity on the Higgs couplings. To this purpose the 
events are subdivided into six different categories in rn™£°. The corresponding numbers are 
reported in Table [V] and will be used in section [V] to extract our bounds on the coefficients 
of the effective operators. 

Let us now discuss the case of a 100 TeV collider. We start by considering a strategy 
similar to the one adopted at 14 TeV, and differing from this one mostly for the numerical 
values of the cuts and selection parameters. The nominal p? threshold on the reconstructed 
jets is increased to pr(j) > 35 GeV to take into account the busier environment of the 
higher energy collision. We define the same pseudorapidity region for the acceptance of jets, 
photons, and leptons although more signal events fall into the high-1 rj region due to the 
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KT [GeV] 

250 - 400 

400 - 550 

550 - 700 

700 - 850 

850 - 1000 

1000 - 

hh 

2.14 

6.34 

2.86 

0.99 

0.33 

0.17 

77 bb 

7.69 

10.1 

3.35 

1.38 

1.18 

0.59 

1133 

0.66 

0.95 

0.31 

0.16 

0.08 

0.045 

tth 

3.33 

4.53 

1.41 

0.41 

0.16 

0.043 

bbh 

0.20 

0.16 

0.03 

0.0054 

0.0022 

0.00054 

Zh 

0.13 

0.19 

0.067 

0.021 

0.009 

0.0009 


TABLE V: Expected number of events after all cuts at \fs = 14TeV in each of the six m T bb ° 
categories considered, assuming an integrated luminosity L = 3000 fb _1 . The last category is 
inclusive. 


fs = 100 TeV 


After selection cuts of Eq. 
After A R cuts of Eq 


24 


After m r ,l co cuts of Eq. 

oo, 77 ^ 

After hadronic W vetoing 


26 


25 


hh 


500 

390 

303 

279 


77 bb 77 jj tth bbh Zh 


22319 2469 

4819 626 

137 18.2 

135 15.5 


2694 

836 

303 

179 


52 102 

17 86 

6.2 3.2 

5.9 3.2 


TABLE VI: Cut flow for the SM signal and the various backgrounds at y/s = 100 TeV. The 
values correspond to the number of events, assuming an integrated luminosity L = 3000 fb _ . All 
reconstruction efficiencies and branching ratios to the final state under study are included. 


boost along the beam axes (this point will be discussed in detail in section IVA). Similarly, 


in our initial selection of events we increase the px cuts on the 6 -tagged jets and photons 
while keeping the mass windows to be the same: 


(26) 


Pt> (' b ), p T > ( 7 ) > 60 GeV , p T< ( b ), p T> ( 7 ) > 40 GeV , 

60 < rri b l co < 200 GeV, 60 < m™° < 200 GeV. 

After the above cuts we find angular distributions similar to those shown in Fig. [7J and we 


thus apply the same cuts as in Eq. (24) and Eq. (25). The cut flow at a/s = 100 TeV is 
reported in Table |VT| Differently from the 14 TeV case, at this level the dominant background 
is tth , rather than 6677 , due to its steeper increase with the collider energy. Imposing a veto 
on extra hadronic activity is beneficial to reduce this background without strongly affecting 
the signal, as illustrated by the distributions of Fig. [9j We find that a veto on hadronic W’s 
is most efficient to maximize the signal significance, and thus apply it. The final number of 


signal and background events after the W veto is reported in the last row of Table VI, while 


a fit to the total signal rate after all cuts is given in Table IV As a final step, we subdivide 


events into six rri T bb ° categories and report the corresponding numbers in Table VII 
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FIG. 9: Multiplicity of jets (left plot) and of reconstructed hadronic W bosons (right plot) after all 
cuts at \/s = 100 TeV. The solid, dotted and dot-dashed curves denote respectively the SM signal, 
77 bb and tth. 


m T h T [GeV] 

250 

-400 

400 

- 550 

550 

- 700 

700 

- 850 

850 

- 1000 

1000 - 

hh 

27.1 

(29.4) 

116.3 

(125.8) 

74.5 

(81.0) 

35.3 

(38.5) 

15.9 

(17.6) 

9.7 (10.8) 

77 bb 

56.1 

(57.5) 

50.5 

(51.3) 

13.1 

(13.1) 

8.37 

(8.38) 

3.63 

(3.90) 

2.79 (3.07) 

1133 

4.23 

(4.96) 

6.14 

(7.38) 

2.54 

(2.87) 

1.34 

(1.46) 

0.46 

(0.63) 

0.78 (0.90) 

tth 

53.7 

(91.3) 

77.3 

(130.3) 

30.1 

(51.7) 

11.5 

(18.7) 

4.37 

(6.99) 

2.49 (3.58) 

bbh 

2.41 

(2.47) 

2.65 

(2.74) 

0.61 

(0.649) 

0.18 

(0.197) 

0.062 

(0.065) 

0.021 (0.026) 

Zh 

0.70 

(0.70) 

1.31 

(1.32) 

0.67 

(0.674) 

0.34 

(0.353) 

0.10 

( 0 . 10 ) 

0.031 (0.031) 


TABLE VII: Expected number of events after all cuts at \fs = 100 TeV in each of the six 
categories considered, assuming an integrated luminosity L = 3000 fb -1 . The last category is 
inclusive. The numbers in the parenthesis are obtained by removing the veto on hadronic IT’s. 

A. Recovering the boosted topologies 

As we anticipated in the previous discussion, a possible issue with our analysis strategy 
is the loss of sensitivity for the kinematical configurations containing boosted Higgses. Al¬ 
though this effect is not likely to be relevant for the 14 TeV LHC, it can potentially affect 
the gg —> hh searches at a future higher-energy collider. In this section we present a first 
estimate of the improvement that can be achieved in the reconstruction of highly energetic 
events by the use of jet substructure techniques. Notice that the development of a fully 
optimized analysis will most likely require a hybrid strategy that smoothly interpolates be¬ 
tween a traditional jet analysis and one using boosted techniques [48] . Devising such a 
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FIG. 10: Left: pseudorapidity distribution of the most forward bottom quark or photon in the 
SM gg —> hh —> 6677 process at 14TeV (solid curve) and 100 TeV (dot-dashed curve). The 
vertical solid line corresponds to the acceptance region at the LHC. The dot-dashed vertical line 
denotes the acceptance region which should be adopted at 100 TeV to retain the same fraction of 
events of the 14 TeV case. Right: A R distance between two bottom quarks or photons in the SM 
gg —» hh —> 6677 process at 100 TeV. The solid (dot-dashed) curve corresponds to the events with 
rrihh > 1000 GeV ( m^h < 1000 GeV). Both plots are made at the partonic level. 


strategy is beyond the scope of the present work and we postpone it to a future study. In 
this section we instead adopt a “minimal” approach and use either the standard analysis or 
the jet substructure technique in each bin. 

There are basically two different effects that lead to boosted Higgses, namely the boost of 
the whole hh system along the beam axis and the production of a Higgs with high p Both 
effects will become relevant at a future high-energy collider. At 100 TeV the di-Higgs system 
acquires on average a non-negligible amount of momentum along the beam axis and, as a 
result, the events more frequently leak into the high-1r/ region. The situation is illustrated 


by the left plot of Fig. [lOj which suggests that a larger pseudorapidity range for object 
reconstruction and flavor tagging could improve the sensitivity to the gg —> hh process. For 
instance, at the 14 TeV LHC roughly 13% of the partonic signal lays outside the acceptance 


region \g\ < 2.5 (solid vertical line in the left panel of Fig. 10). If one assumes the same 
acceptance region, this fraction increases to ~ 30% at a 100 TeV collider. In order to obtain 
the same acceptance efficiency of the 14 TeV LHC, the coverage region should be extended to 
|?y| < 3.3 at the higher-energy collider (vertical dot-dashed line in the left panel of Fig. 10). 17 


17 After imposing the cuts in Eqs. (23) and (241 on partonic signal events, these fractions reduce to ~ 6.5% 
for the 14 TeV LHC and ~ 24% for the 100 TeV case. To recover the 14 TeV acceptance efficiency at the 
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A further problem with boosted events is that our analysis relying on a dijet-style search 
fails when the two 6 -jets are not well separated. A similar argument can be applied to the 
di-photon system. The boosted Higgses are more likely to be produced in the high invariant 


mass tail of the m^h distribution. This effect is illustrated in the right panel of Fig. 10 


where we plot the A R separation between two 6 ’s or photons at the partonic level for the 
two categories of events with rrihh > 1000 GeV and rrihh < 1000 GeV. Most events belonging 
to the lower rrihh category have a clear two-prong topology, whereas a significant amount of 
events in the higher m hh category fail to be resolved as two well separated partons. 

For our present purposes it will be sufficient to perform a simple analysis on the 66 
subsystem and compare the performance of our traditional jet-based analysis with the sub¬ 
structure method. For the 77 subsystem we consider two possibilities: a reduced isolation 
cone size (Ri SO = 0 . 2 ) and an ad-hoc mutual photon isolation (a somewhat similar prescrip¬ 
tion was proposed in |49] for a boosted di-tau system). In the latter, photon isolation is 
imposed neglecting the other photon in the 77 system and an event is retained as long as 
Ai 2 ( 7 , 7 ) > 0.2. This procedure is quite similar to our traditional jet-based analysis up 
to the slight modification of the photon isolation criterion. We do not modify instead the 
isolation criterion on leptons. 18 

After this step, the events are clustered into R = 1.5 “fat-jets” with the Cam¬ 
bridge/Aachen (C/A) jet algorithm [ 501 15T] . We iteratively look over those fat-jets and 
apply the BDRS sub jet-finding technique of Ref. [52] (with the same declustering parame¬ 
ters as in [52]). To ensure enough boost we only look into fat-jets with pr(j) > 150 GeV. The 
declustering algorithm stops when it successfully identifies three subjets (the third subjet 
takes into account the leading gluon emission from either bottom quark line). The fat-jet 
is identified as a Higgs-jet if the invariant mass of the three subjets falls into the Higgs 


mass window mj| co = [105,145] GeV (see Eq. (25)). If multiple candidates exist, we pick 
the one whose invariant mass is closest to the Higgs mass. Two 6 -taggings are performed 
at subjet-level assuming a 70% 6 -tagging efficiency. The cuts on the angular separations in 


Eq. (24) are applied to the subjets and the photons. Finally, we do not apply any veto on 
extra hadronic activity or on hadronic W bosons. 

The performances of the substructure analysis on signal events are summarized in Ta¬ 


ble |VIII| and compared to the traditional jet-based analysis. We consider three possible 
scenarios which differ by the treatment of the 77 system: “Substructure I” uses a stan¬ 
dard photon isolation criterion with R iso = 0.4, whereas in “Substructure II” the cone size 


higher-energy collider the pseudorapidity region should be extended to |? 7 | < 3.6. 

18 While naively one would expect that leptons from the boosted tops in the resonant tth background fail 
more frequently the isolation requirement, we find that the efficiency drop is not significant for the mhh 
range of interest. 
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Signal event rate at yfs = 14 TeV, assuming L = 3 ab 1 


m T ti° [GeV] 

250-400 

400-550 

550-700 

700-850 

850-1000 

1000- 

Trad. anal. 

2.14 

6.34 

2.86 

0.99 

0.33 

0.17 

Substr. I 

0.21 

3.53 

2.31 

0.81 

0.31 

0.22 

Substr. II 

0.25 

3.91 

2.52 

0.88 

0.34 

0.24 

Substr. Ill 

0.21 

3.53 

2.31 

0.81 

0.32 

0.23 


Signal event rate at \/s = 100 TeV, assuming L = 3 ab 

-1 


mlT [GeV] 

250-400 

400-550 

550-700 

700-850 

850-1000 

1000-1150 

1150-1300 

1300-1450 

1450-1600 

1600- 

Trad. anal. 

29.3 

125.6 

80.8 

38.4 

17.6 

7.04 

2.27 

0.68 

0.52 

0.23 

Substr. I 

4.1 

70.6 

64.7 

32.1 

16.6 

9.64 

4.02 

1.90 

1.21 

0.99 

Substr. II 

5.39 

89.8 

78.1 

37.5 

19.3 

10.5 

5.0 

2.81 

2.05 

4.7 

Substr. Ill 

4.1 

70.6 

64.7 

32.5 

17.0 

9.71 

4.40 

2.50 

1.82 

4.02 


TABLE VIII: Comparison between the traditional jet-based analysis (second line) and the substruc¬ 
ture analysis (third, fourth and fifth lines) on SM signal events. The three scenarios “Substructure 
I, II, III” differ by the treatment of 77 system and are defined in the text. 


is reduced to Ri SO = 0.2. “Substructure III” corresponds instead to the mutual isolation 
criterion. Not surprisingly, the substructure analysis becomes more efficient compared to 
the traditional one for m r ^° > 1000 GeV; this matches our naive expectation based on the 
A R(b,b) and Ai 2 ( 7 , 7 ) distributions of Fig. 10 For < 1000 GeV, on the other hand, 
only part of the events are correctly reconstructed using the substructure technique (due 
to the finite fat-jet size R = 1.5) and the traditional analysis is more efficient. This is also 
clearly illustrated by the plots of Fig. 11, which show the mjy°° distributions for signal events 
with both strategies. These results suggest that the use of jet substructure in the 6677 final 
state is crucial only at a 100 TeV collider, while it has a minor impact at the 14 TeV LHC, 
where large values of rrihh are suppressed by the fast drop of the gluon luminosity. The 
relevance of jet substructure can be determined in a more quantitative way by comparing 
the reach in at 14 TeV and 100 TeV. At the partonic level, Fig. [ 6 ] shows that the same 
number of events one has at 14 TeV after imposing the cut rrihh > 850 GeV is obtained 
at 100 TeV for rrihh > 2570 GeV. At the hadronic level, the increase in the reach is much 
more modest than the naive expectation when adopting a traditional jet-based analysis. For 
instance, the signal rate at 14 TeV is N events ~ 0.5 for m^° > 850 GeV. The same number 
of events is obtained at 100 TeV for m^° > 1540 GeV using a traditional analysis, and 
> 2560 GeV using jet substructure with the reduced photon isolation cone R iso = 0.2. 
Hence, a jet substructure analysis (along with the appropriate modification of the photon 
isolation criteria) is essential to get close to the naive expectation. 

Focusing on the 100 TeV case, we show in Table IX the number of signal and background 


events expected with the substructure analysis and a reduced photon isolation cone size 














Signal (SM), 14TeV | 


| Signal (SM), 100 TeV | 




FIG. 11: Number of signal events as a function of mj ,^ 0 after all cuts at 14 TeV (left plot) and 
100 TeV (right plot). The solid line is obtained by adopting the substructure technique with 
reduced photon isolation cone size (scenario “Substructure II” described in the text), whereas the 
dot-dashed line corresponds to a traditional jet-based analysis. 


KT [GeV] 

850 - 1000 

1000 - 1200 

1200 - 1400 

1400 - 1600 

1600 - 1800 

1800- 

hh 

19.3 

12.5 

4.86 

3.03 

1.75 

2.96 

77 bb 

3.92 

2.35 

1.18 

0.59 

0.29 

0.098 

tth 

8.63 

4.91 

2.32 

1.16 

0.51 

1.1 


TABLE IX: Expected numbers of signal and background events at y/s = 100 TeV after the jet 
substructure analysis (scenario “Substructure II” described in the text with photon isolation cone 
Riso = 0.2) assuming an integrated luminosity L = 3000 fb -1 . The last category is inclusive. 


(scenario “Substructure II") in six categories with high m™£°. We include only the dominant 
6677 and tth backgrounds for simplicity. 19 A further reduction of the tth background could 
in principle be possible by iteratively looking for boosted top- and W -jets through dedicated 
tagging algorithms, and by discarding events where such objects are reconstructed. This is 
however beyond the scope of the present work. On the other hand, we have checked that a 
simple veto on hadronic W's reconstructed from pairs of jets, as applied in the traditional 


19 We have checked that Zh and bbh are negligible, as the corresponding number of events in each category 
is always smaller than 0.1. We expect that also the background 77 jj can be reduced to a subdominant 
level as long as the efficiency for making two &-tags at the subjet level is sufficiently high. The 77 bb 
background was newly generated for the substructure analysis with relaxed generation cuts A R(b, b) > 0 . 1 , 
Ai?(y, 7 ) > 0.15; these are less stringent than those specified in Appendix 
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analysis, is not efficient on events with such high in T ^°. The results in Table IX will be used 
in Section [V] to assess the impact of jet substructure in setting bounds on the coefficients of 
the effective operators. 


V. SENSITIVITY ON THE EFT COEFFICIENTS 

The results of the previous section have been processed through a Bayesian statistical 
analysis to extract the sensitivity on the coefficients of the Higgs effective Lagrangian. In 
each case, a probability distribution of the relevant parameters is obtained by constructing 
a likelihood function from the signal and background number of events in the six m^h cate¬ 
gories, and marginalizing over (or fixing) the remaining parameters. The injected signal is 
the SM one for all the results presented in this section. A flat prior is assumed for the coef¬ 
ficients of the Higgs effective Lagrangian, except when they are constrained by single-Higgs 
measurements. In the latter case, we use the ATLAS projections for the high-luminosity 
LHC of Ref. [53] to construct a likelihood function and use it as a prior. Reference [53] 
reports the estimated precision on various Higgs decay channels expected for a/ s = 14 TeV 
with integrated luminosities L = 300 ftY 1 and L = 3ab _1 . We will further study the case of 
a future proton-proton circular collider operating at a/s = 100 TeV with 3 ab _1 of integrated 
luminosity. When single-Higgs measurements are required to perform marginalization in 
this latter scenario, we use the ATLAS projections for L = 3ab _1 at a/s = 14 TeV, lacking 
a specific estimate of the precision reachable at a 100 TeV machine. We will thus consider 
the following three benchmark scenarios: 

LHC 14 HL-LHC FCCioo 

a/s 14 TeV 14 TeV 100 TeV 

Luminosity L = 300fb~ 1 L = 3ab -1 L = 3ab -1 

For simplicity, we do not include in our analysis the theoretical uncertainty on the predic¬ 
tion of the signal cross section nor the systematic uncertainties which will characterize the 
extraction of the background from data in a realistic experimental analysis. 20 

As a first result, we derive the precision that can be obtained on the signal strength 
multiplier /i = cr/asM hi the three benchmark scenarios described above. No marginalization 
is performed in this simple case. We obtain the following 68 % intervals on /i\ 

20 The theoretical error was estimated by the authors of Ref. [7] to amount to a ~ 18% (12%) from scale 
variation, ~ 10% (10%) from the use of EFT in the calculation of the NLO k-factor, and ~ 7% (6%) 
from PDFs at the LHC 14 (FCCioo)- See Ref. [53] for a proposal on how to reduce the theoretical and 
experimental uncertainties by considering the ratio of double and single cross sections. 
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FIG. 12: 68 % probability contours in the planes (c 2 t,c 3 ) (left plot) and (c 2 t,c 2s ) (right plot). The 
different curves refer to the three benchmark scenarios: LHC 14 (black continuous line); HL-LHC 
(black dashed line); FCC 100 (red dotted line). 



LHCm 

HL-LHC 

FCC 100 

68 % interval on /i 

[-0.41,3.0] 

[0.50,1.6] 

[0.92,1.1] 


Very similar results are obtained from an inclusive analysis without rrihh categories, as one 
naively expects since /i accounts only for an overall rescaling of the total cross section of the 
signal. 

Turning to the determination of the coefficients of the effective Lagrangian, we first 
consider the case of the non-linear parametrization of Eq. Q. Figure 12 shows the 68 % 
probability contours in the planes (c 2 t ,c 3 ) and (c 2 t,c 2s ) obtained through the procedure 
described above. In this case the marginalization is performed over two parameters: c t , with 
a prior obtained from single-Higgs measurements, and the branching ratio for the decay 
hh —y bb'y'f. For the latter parameter we assume a Gaussian distribution around the SM 
value with standard deviation equal to 0.15 and 0.10 respectively for L = 300fb _1 and 
L = 3ab -1 (both at 14TeV and 100 TeV). For simplicity, the remaining couplings are set 


to their SM values: c g = c 2g = 0 in the plot on the left of Fig. 12 c 3 = 1 and c g = 0 in the 
plot on the right. We have checked that performing an additional marginalization over c g 
slightly decreases the precision on the measured couplings, without changing the shape of 
the contours. 

We fold that the couplings c 3 and c 2t are strongly anti-correlated, and the precision 
expected on c 2 * is much higher than the one on the Higgs trilinear coupling c 3 . This is 
in agreement with previous studies, which pointed out the strong sensitivity of the dou- 
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FIG. 13: 68 % probability contours for the FCCioo with a traditional analysis (black continuous 
line) and one using jet substructure (red dotted line), see text. Left plot: plane ( 024 , 03 ); Right 
plot: plane (c 2 t ,c 29 ). 


ble Higgs cross section on the tthh quartic coupling, see [ 12 ] and references therein. The 
coupling C 2 g is determined even more accurately, although its naive estimate is suppressed, 
compared to that of C 3 and c 2t , by two powers of a weak spurion if the Higgs is a pseudo 
NG boson, see Eqs. (§, Although the Higgs trilinear coupling c 3 is the less accurately 
measured parameter, its precision can be highly increased at a 100 TeV collider, as shown 
by the left plot of Fig. [l2j This is mainly due to the higher number of signal events which 
can be produced at this machine. Using boosted jet techniques does not seem to improve 
significantly the accuracy on C 3 and c 2t , while it increases dramatically that on c 2ff . This is 
expected, since jet substructure techniques are most relevant to reconstruct highly boosted 
events with large rrihh, which are especially important to determine c 29 . The situation is 
illustrated by Fig. 13, where the red dotted curve is obtained by including the first 5 cat¬ 
egories (250 GeV < rrihh < 1000 GeV) of the traditional jet-based analysis and the last 5 


categories (rrihh > 1000 GeV) of the “Substructure 11” analysis, see Tables VII and IX 


In order to break the degeneracy among the various parameters and extract them 
precisely, it is crucial to make use of the information on rrihh■ We find that an inclusive 
analysis, where events are not classified in the six rrihh categories of Tables [V] and VII 


is much less powerful in constraining the Higgs couplings, especially in the case of a 
100TeV collider. This is illustrated by the plots of Fig. 14 in the plane (c 2 t ,c 3 ). It is 
evident how at 100 TeV only an exclusive analysis can break the degeneracy between 


c 2 t and C 3 . As expected from the discussion of Section IIC, categories with larger rrihh 
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C2t 


C2t 


FIG. 14: 68 % probability contours in the plane (c 2 t, C 3 ) with the exclusive analysis (black continuous 
line) and an inclusive one without m^h categories (red dotted line), see text. Left: plot for the 
HL-LHC; Right: plot for the FCCioo- 


most strongly constrain c 2t and c 2g , while C 3 is mainly determined by events at threshold, 
see Fig. [l5j Given the relevance of the categories with large rrihh in determining the 
Higgs couplings, it is important to give an assessment on the validity of the effective field 


theory approach in our analysis. Following the discussion of section IIB , we evaluate 
the minimum coupling strength g m i n ~ VS(E/v) by estimating the maximal invariant 
mass E of the events which lead to a precision 5 on c 2ff (i.e. the parameter controlling 
the term which grows quadratically with the energy in the scattering amplitude). To 
this aim we re-derive the 68 % probability contours in the plane (c 2 t, c 2s ) by removing 
one or more of the categories with large rrihh, as a way to estimate their impact on the 
determination of c 29 . We show such plots in Fig. [l 6 ]for the HL-LHC and the FCCioo- We find 


LHC 


14 


HL-LHC 


FCC 


100 


FCCioo (boosted) 


5 ~ 0.15 
E ~ 1.0 TeV 

9min 1-6 


5 ~ 0.1 

E ~ 1.0 TeV 

9min 1-3 


5 ~ 0.05 
E ~ 1.0 TeV 
9min 0.9 


5 ~ 0.016 
E ~ 1.8 TeV 
Qmin 0.9 


where the last column refers to the analysis including jet substructure at the FCCioo (first 
5 categories of traditional analysis + last 5 categories of “Substructure II” analysis). As 
discussed in Section IIB , the value of g m i n determines the extension of the region which 
can be probed under the validity of the effective field theory, see Fig. [T[ In particular, 
gmin < Ut on ly in the case of a 100 TeV collider, which suggests that an analysis including 
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FIG. 15: Breakdown of rrihh categories in the plane (c 2 t,cs) for the HL-LHC (left plot), and in 
the plane ( C 2 t,C 2 g ) for the FCCioo (right plot). The various curves indicate the 68% probability 
contours for the following pairs of categories of Tables |V| and VII 1 and 2 (dotted blue line); 3 and 4 
(dot-dashed brown line); 5 and 6 (dashed red line). 
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FIG. 16: 68% probability contours in the plane (c 2 t, C 2 g ) for the HL-LHC (left plot) and the FCCioo 
(right plot). The different curves have been obtained by removing the following rrihh categories of 
Tables |V| and VII from the fit: 6 (dashed red line); 6 and 5 (dot-dashed brown line); 6, 5 and 4 
(dot blue line). The continuous black contour is obtained by including all the categories in the fit. 


34 











only dimension -6 operators at the LHC (even at its high-luminosity upgrade) is not sensitive 
to the case of a pseudo NGB Higgs with fully or partially composite Ir, although it can probe 
the case of a generic composite Higgs. 

We now turn to the effective Lagrangian for a Higgs doublet, Eq. (J2]), and show the 
constraints on its coefficients. The production cross section, in this case, depends on ch , 
c u , c g and Cq, while the branching ratio for hh —* 6677 depends on ch , c u , c g as well as on 
c ( i and c 7 . The latter two coefficients parametrize respectively the modihcation of the down 
quark Yukawa couplings, in particular that of the bottom quark, and the contact interaction 
between the Higgs boson and two photons. Their precise definition can be found in Ref. na 
and it is analogous to that of c u and c g in Eq. (J 2 J. We will set c 7 = 0 for simplicity in the 
following. 

We start considering the constraints on c u and c g that follow from both single and dou¬ 
ble Higgs production. The 


probability contours are shown in Fig. [17| for the three 
benchmark scenarios considered in this section. The solid blue curve refers to double Higgs 
production in the 6677 final state (our analysis), while the dashed and dotted blue curves 
correspond to the constraint from respectively all single-Higgs processes except tth, and 
tth alone 21 (the ATLAS projections from Ref. [53]). They have been obtained by fixing 
ch = c d = 0 and marginalizing over c 6 with flat prior (in the case of double Higgs pro¬ 
duction). The two dashed ellipses in the plots in the upper row correspond to the two 
degenerate solutions that follow from the interference between c g and c u in the gg —h cross 
section. The filled areas denote instead the 68 % probability regions obtained by combining 
the following processes: all single Higgs processes including tth (orange region); all single 
Higgs processes plus double Higgs production (dark blue region); all single Higgs processes 
except tth plus double Higgs production (light blue area in the bottom left plot of Fig. [Tf| ). 
They have been obtained by marginalizing over c#, c d and Cq. One can see that double 
Higgs production breaks the degeneracy on c g that remains even after including tth among 
single-Higgs processes, ft also helps increasing the precision on both c u and c g compared 
to that following from single Higgs processes alone, especially at a 100 TeV collider. I 11 this 
respect, it is interesting to notice that c u and c g are strongly correlated when considering 
single-Higgs measurements without tth, and that the extension of the dashed ellipses is much 
larger than the error obtained on each individual parameter by setting the other to its SM 
value (in fact, the ellipses become infinite bands if one lets c 7 vary, as a consequence of 


21 Notice that the dotted curves are not vertical in the plane (c u ,c g ) due to the dependence of the total 
Higgs width on c g . A further dependence on c g would follow from the contribution of diagrams with an 
insertion of O g to the tih cross section; this effect has not been included in our fit. We thank C. Grojean 
for drawing our attention on this point. 
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FIG. 17: 68 % probability contours in the plane (c u ,c g ) for the three benchmark scenarios: LHC 14 
(upper left plot); HL-LHC (upper right and lower left plots); FCC 100 (lower right plot). Different 
curves correspond to: double Higgs production in the 6677 final state (solid blue line); all single- 
Higgs processes except tth (dashed blue line); tth alone (dotted blue line). Filled areas correspond 
to: all single Higgs processes including tth (orange region); all single Higgs processes plus dou¬ 
ble Higgs production (dark blue region); all single Higgs processes except tth plus double Higgs 
production (light blue area in the lower left plot). 


the degeneracy in the h — > 77 decay rate) 22 . The process tth mainly constrains c u and is 
crucial to reduce the overall experimental uncertainty. In comparison, gg —> hh —> 6677 is 
less powerful but still competitive in constraining c u . For example, by removing tth from the 

22 By fixing all the parameters to their SM value but one, the likelihood obtained from Ref. [53] is approx¬ 
imately Gaussian. By including all single Higgs processes except tth we obtain the following standard 
deviations: cr(c u ) = 0.06, a(c g ) = 0.005, ct(ch) = 0.08 for L = 300fb _1 ; cf(c u ) = 0.05, a(c g ) = 0.004, 
cj(ch) = 0.05 for L = 3ab _1 . From tth alone we find: a(c u ) = 0.2 for L = 300 fb _1 ; a(c u ) = 0.08 for 
L = 3ab _1 . 
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FIG. 18 : Left plot: 68% probability contours in the plane {c u ,cq)\ Right plot: probability distri¬ 
butions as functions of cq. The different curves refer to the following benchmark scenarios: LHC14 
(black continuous line); HL-LHC (black dashed line); FCC100 (red dotted line). 


combination of single and double Higgs at the HL-LHC, one obtains the light blue region 
(instead of the orange one) which is sensibly smaller than the dashed contour. One must 
also notice that the tth dotted curve in Fig. 17 corresponds to the combination of several 
decay channels (all those studied by Ref. [53]: h —» 77 ,ZZ,nn), while the continuous blue 
curve denotes double Higgs production in the 6677 final state alone, i.e. the one studied in 
this work. Combining with other final states (e.g. bbrr and bbWW) will make double Higgs 
production even more competitive for the determination of c u , as well as of c g . 

The plot on the left of Fig. 18 shows the 68 % probability contours in the plane (c u ,c 6 ), 


obtained by marginalizing over ch , Q and c g . As already suggested by Fig. 12 , the precision 
on cq (i.e. the parameter which controls the modification of the Higgs trilinear coupling) is 
much smaller than that on c u . At i/s = 14TeV, in particular, the constraint is on values of 
Cq larger than 1. 23 Further marginalizing over c u leads to the probability functions shown in 


the right plot of Fig. 18 for the three benchmark scenarios. Notice that even at the HL-LHC 
the likelihood is far from being Gaussian, and a second maximum is present at Cq ~ 4.5. We 
find the following 68 % probability intervals on Cq: 


LHCm HL-LHC FCC 100 

[-1.2,6.!] [-1.0,1.8] U [3.5, 5.1] [-0.33,0.29] 


23 See the discussion of Section 


IIA 


on the validity of the effective field description in this case. 
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FIG. 19: Probability distributions as functions of cq at the HL-LHC. The black continuous curve 
refers to the standard exclusive analysis discussed in section IY The dotted red curve refers to an 
inclusive analysis in the left plot, and to an analysis without marginalization on ch, Cd, c g and c u 
in the right plot. See text. 


It is interesting to compare the probability functions of Fig. 18 with those obtained with 
an inclusive analysis without rrihh categories. For the HL-LHC for example, an inclusive 


analysis gives the function shown in the left plot of Fig. 19 (dotted red curve). The two 
solutions are in this case almost degenerate, and in fact the second (non-SM) maximum 
is slightly higher than the one at Cq = 0 as an effect of the marginalization on c u 24 . A 
fully exclusive analysis thus removes the degeneracy and helps reduce the significance of the 
unphysical solution. 

Another interesting question concerns the relevance of the marginalization over c#, c<f, 


c g and c u in determining the Higgs trilinear coupling. The right plot of Fig. 19 shows the 
probability function which is obtained by switching off the marginalization in the HL-LHC 
scenario (dotted red curve). The effect of marginalization is that of flattening the SM 
solution, although in a marginal way at 14 TeV. The effect is on the other hand much more 
important at 100 TeV. On a more quantitative side, the 68 % probability intervals which 
follow without marginalization are: [—1.4, 5.8] (LHC 14 ), [—1.0,1.6] U [3.8, 5.1] (HL-LHC), 
and [—0.18,0.18] (FCCioo)- It is clear that a precise measurement of the Higgs trilinear 
coupling, as it follows at 100 TeV from the 6677 channel or as it might be obtained at 
the high-luminosity LHC from a combination of several decay modes, relics of the accurate 
extraction of the other couplings. In particular, a significant uncertainty on cq could follow 


24 This “unwanted” feature can be avoided by profiling, instead of marginalizing, over c u . The second 
maximum in this case would be lower than the one at Cq = 0, since the highest peak of the 2-dimensional 
likelihood is indeed the one at the SM point. 
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FIG. 20: Width of the 68% probability interval on c6 (orange region) as function of the uncertainty 
on c u for the FCCioo- See text. 


from a poor determination of c u . This is illustrated in Fig. 20 for the 100 TeV scenario. 
The plot shows the extension of the 68% probability interval on cq as a function of cr(c u ), 
where c u is marginalized with a Gaussian distribution with standard deviation a(c u ). The 
uncertainty on Cq grows with a(c u ) for small values of this latter parameter, while it becomes 
constant for a(c u ) > 0.2. In this limit double Higgs production alone determines c u with a 
smaller error, independently of single-Higgs measurements. 

The impact of the uncertainty on the top Yukawa coupling in a precise determination of 
the Higgs trilinear was already discussed in Ref. [53], although in a scenario where the only 
modifications with respect to the SM are in the values of these two couplings. Such scenario 
of New Physics, however, is a disfavored one. If the Higgs boson is part of a weak doublet 
and the new states are heavy, one gets the effective Lagrangian (|2]) at low energy, where the 
same operator O u that gives a modification to the top Yukawa coupling also generates the 
quartic tthh interaction, whose contribution must thus be included. On the other hand, it is 
possible to get a modification of the top Yukawa without having a direct tthh interaction if 
the new states are light, as for example in two-Higgs doublet models. In that case, however, 
the new states will give an extra contribution to the production cross section, which is in 
fact required to match that from the tthh vertex in the decoupling limit. One can envisage 
a situation where the top Yukawa coupling c t is modified by new heavy physics while c 2 t = 0 
if the Higgs does not belong to a weak doublet, like for example in the case of a dilaton. 
This class of models is however disfavored by current data. These considerations suggest 
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that any New Physics which modifies the top Yukawa coupling is likely to have an additional 
impact on the double Higgs cross section. More specifically, if the new states are heavy and 
the Higgs is part of a doublet, then any uncertainty on the value of the top Yukawa coupling 
will reflect on a contribution from the tthh vertex which must be properly included. Of 
course one can restrict oneself to the SM case and ignore the tthh anomalous coupling; in 
that case however, the top Yukawa coupling can be more precisely determined from the top 
quark mass, and its uncertainty is sufficiently small to have a negligible impact on double 
Higgs production. 

VI. COMPARISON WITH PREVIOUS WORKS 

The results of the previous section suggest that an analysis of double Higgs production 
in the 6677 final state can determine the Higgs trilincar coupling with a ~ 30% precision 
at the FCCioo- 25 At the HL-LHC, on the other hand, the estimated uncertainty is at 
the level of 200 %, and in fact another solution of the fit is present for a trilinear coupling 
approximately equal to 5.5 times its SM value. These results are more pessimistic than other 
estimates that appeared in the literature, and it is thus useful to compare our analysis with 
previous works in order to highlight the different strategies and assumptions. 

The first detailed analysis of the 6677 final state was performed by Baur, Plehn and 
Rainwater in Ref. [3]. They retained events with one or more 6 -tags at the LHC, while two 6 - 
tags were required at its high-luminosity phase. An analysis with only one 6 -tag seems much 
more involved due to the larger background and the systematic error from combinatorics. 
We will thus compare the results for the HL-LHC obtained with two 6 -tags. A first difference 
with our analysis is in the value of the 6 -tagging efficiency and the corresponding light jet mis- 
tag rate, since Ref. [4] assumes = 0.5 and = 1/23, to be compared with e& = 0.7 and 
ej^.b = 0.01 used in the present work (for both the LHC and its high-luminosity upgrade). 
While the choice of Ref. [4J was based on early Montecarlo simulations, more recent results 
based on LHC data have shown that larger 6 -tagging efficiencies can be achieved while 
maintaining acceptable levels of rejection rates [4f)fj42] . As a consequence of the larger value 
of €j^b, the background jj 77 is found to be the largest one in Ref. [4], differently from 
our estimates. A second difference concerns the estimate of the background 6677 , which we 
found to be enhanced at NLO by a large k-factor k bbll ~ 2. Reference [4] instead computes 
all the backgrounds at LO and rescales their cross sections by an ad-hoc common factor 1.3 
to account for NLO effects, thus underestimating 6677 . Finally, a narrower mass window 

25 Here and in the following, by sensitivity/precision on some Higgs coupling or coefficient of the effective 
Lagrangian we mean the 68% error on its measured value for injected SM signal. 
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for the photon pair, m r ff° = ± 2.3 GeV, is assumed in the study of Ref. [4J. This leads to 

a reduction by a factor ~ 2 compared to our estimates of all the non-resonant backgrounds, 
since their distribution is nearly flat around the Higgs mass. We have chosen instead 
to make a more conservative cut on m T ff° and to leave the issue of optimizing the width of 
the mass window to a fully-realistic experimental analysis. As a consequence of the above 
different assumptions, the total background in Ref. [4] is rather smaller than in our analysis, 
which led the authors to a more optimistic estimate of the precision on the trilinear coupling. 

A more recent analysis of the 6677 channel has been carried out in Refs. d ei mg. 
Reference [7] by Baglio et al. estimates 47 SM signal events after all cuts with L = 3ab^ x 
at the high-luminosity LHC, which should be compared with our ~ 13 events. Although 
our smaller rate can be in part explained by the inclusion of the photon efficiency factor 
0.8 2 = 0.64 (not included in Ref. [7]), we were not able to fully identify the origin of such 
difference. Furthermore, the 6677 continuum is found to be negligible in Ref. [7] due to 
a limited MC statistics, thus leading to an underestimation of the background rate. The 
comparison with Ref. [TOj is instead difficult, since the number of signal and background 
events is not reported. The same kinematic cuts of Ref. [7J are however applied, and it is 
said that an agreement is found with the efficiencies there reported. This suggests that the 
same underestimation of the background also plagues Ref. m. 

Finally, we compare with the analysis by Yao in Ref. [ 8 ], upon which the Snowmass 
study fS! is based. This is a sophisticated study including a full-fledged simulation of 
showering, hadronization and detector effects. The estimated number of SM signal and 
background events at the LHC after all cuts is respectively 16.6 and 53.4 (40.1 of which 
from 6677 ) for L = 3ab _1 , with a corresponding statistical significance S/\fB = 2.3. This 
is compatible with the numbers reported in Table III which give S/\fB = 2.1. Our analysis 
thus agrees with the signal and background rates estimated by Yao. Based on these numbers, 
this latter calculates a precision on the trilinear coupling equal to 50% and 8 % respectively for 
the HL-LHC and the FCCioo- This is much more optimistic than our corresponding estimates 
200 % and 30% reported in Section |VJ The discrepancy follows because the uncertainty on 
the trilinear is derived in Ref. ( 8 j from the dependence of the signal cross section on this 
parameter before cuts. In particular, a linear approximation is used around the SM point 
with a slope d(cr/ gsm) / dc 3 ~ —0.8 taken from Ref. [7j. However, the slope decreases 
substantially (in absolute value) after imposing all cuts: at the end of our analysis we find 
d(a/a SM )/dc 3 ~ —0.58 and —0.50 respectively for the HL-LHC and the FCCioo- 26 By 
using these numbers, the uncertainty on the trilinear coupling that follows from the rates 


26 In fact, the dependence on the trilinear coupling also depends on the rrihh category considered, the 
categories with lowest rrihh having the largest absolute slope, as expected. 
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of Ref. [ 8 ] is equal to 90% and 16% respectively for the HL-LHC and the FCCqoo- 2 ' The 
difference between these improved values and our results is finally due to the fact that the 
linear approximation is not accurate for y/s = 14TeV (a second solution exists), and to 
the marginalization on c u , ch, Cd, c g (not included in Ref. [ 8 ]) which has a large effect at 
yfs = 100 TeV as explained in Section [vj 

Besides the above theoretical studies, the ATLAS collaboration has recently completed 
a Montecarlo analysis of the 6677 decay mode for the HL-LHC [24J. The simulation of the 
6677 background is done by matching up to 1 extra jet at the matrix-element level, thus 
properly including the bulk of the NLO corrections. The reported signal and background 
rates are consistent with those obtained in our analysis, although backgrounds with fake 
photons ( btrfj) and fake 6 -jets from charm quarks (CC 77 ) are found to be sizable, differently 
from what assumed in this work. It is interesting to see if further experimental analyses 
will confirm this preliminary study and strategies will be found to reduce the size of these 
reducible backgrounds. 


VII. SUMMARY AND OUTLOOK 

In this work we presented an analysis of double Higgs production via gluon fusion. The 
novelty of our approach with respect to most of the previous studies is the use of an effective 
field theory perspective that allows one to encode all possible effects from heavy New Physics 
into a small set of deformations of the SM Lagrangian. Although the leading contributions to 
observables are in general expected to arise from dimension -6 operators, when these latter are 
suppressed due to selection rules higher-dimensional operators can become important. We 
pointed out that this occurs in double Higgs production, where the dimension -6 operator O g 
violates a Higgs shift symmetry and is thus suppressed if the Higgs is a pseudo Nambu- 
Goldstone boson. I 11 this case dimension -8 operators become relevant in the high invariant- 
mass tail of the kinematic distributions. A careful assessment of the range of validity of the 
effective field theory description is thus required, like the one provided in Section |TT| There 
we analyzed in detail the implications of the standard SILH power counting, which predicts 
0{y 2 /f 2 ) deviations in a large set of observables, including the Higgs trilinear coupling. 
We discussed alternative scenarios, including a Higgs-portal model, that imply a modified 
power counting and lead to large deviations in the trilinear coupling while predicting small 
modifications to the other Higgs couplings. These will be presumably the first scenarios 
to be probed by early measurements of double Higgs production at the LHC, where the 

2l During the completion of this paper, an updated analysis by Yao was presented in a conference talk [SB] , 
where a similar precision is obtained at 100 TeV after taking into account the correct slope. 
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precision on the trilinear coupling is expected to be limited. 

We provided an explicit application of the effective held theory approach by assessing the 
experimental sensitivity on the coefficients of the relevant local operators. For def ini teness 
we focused on the gg —>■ hh —>• 77 bb process, which has been recognized as one of the cleanest 
channels to exploit double Higgs production. Although this process was previously analyzed 
in several works, our study provides some important improvements which led to significantly 
different results (see Section VI for a full comparison). One of the most relevant observables 
that can be extracted is the Higgs trilinear coupling C 3 . We found that the 77 bb channel 
allows the determination of C 3 with a fair precision (~ 30% with L = 3 ab _1 ) only at a future 
100 TeV hadron collider. The prospects for the LHC, instead, are much less optimistic and 
only an 0(1) determination seems possible even with 3 ab -1 of integrated luminosity. This 
result is significantly worse than what usually claimed in the literature. The origin of the 
discrepancy is in large part due to a more careful determination of the background processes. 
In particular we found that the irreducible 6677 background is enhanced by a sizable NLO 
k-factor ( k b i 77 ~ 2 ), not included in most of the previous works, and is thus larger than 
previously thought. We also showed that the determination of C 3 is affected by the choice 
of the statistical treatment. I 11 particular a full fit including possible corrections to all 
relevant dimension -6 operators leads to significant differences with respect to the procedure, 
often used in the literature, where only the Higgs trilinear coupling is allowed to vary. An 
important role is played by the uncertainty on the coupling of the Higgs boson to the top 
quark. The reason is that, in the context of the effective Lagrangian for a Higgs doublet, the 
same operator O u which controls the modification to the top Yukawa coupling also generates 
a quartic tthh anomalous interaction which leads to a new contribution to the double Higgs 
production cross section. This latter effect must be properly included when estimating the 
precision on the Higgs trilinear coupling. We find that an uncertainty of order 10% on 
the measurement of c u can nearly double the error on C 3 at a future 100 TeV collider (see 
Fig. 20). One can in fact turn the argument around and use double Higgs production to 
determine c u . Our results show that the accuracy one can obtain in this way is potentially 
competitive with the one from the tth process (see Fig. 0 . especially if several final states 
are combined. On more general terms, other couplings can be extracted from the gg —y hh 
process. In particular, producing a pair of Higgses gives unique access to the tthh and gghh 
quartic interactions, parametrized respectively by c 2 t and c 2g . These should be regarded as 
independent parameters from the wider perspective of a non-linear effective Lagrangian. We 
found that these couplings can be determined much more accurately that C 3 , as shown in 
Fig. [I2j in agreement with previous works. 

Another important point of our study is the effectiveness of an exclusive analysis exploit¬ 
ing the rrihh differential distribution. This procedure is particularly useful to disentangle 
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the contributions arising from different effective operators since each of them leads to very 
specific deformations of the rrihh distribution. This analysis strategy can be relevant for 
the LHC mainly in its high-luminosity phase, and can have a dramatic impact at a future 


100 TeV collider (see Figs. 14 and 15). The other kinematic variable characterizing double 
Higgs events, namely the angular separation between the two Higgses, was instead found to 
have a marginal role in our analysis. The expected distribution, indeed, is almost flat in the 
SM and nearly unchanged by New Physics effects. A possible exception to this result are the 
deformations due to dimension -8 operators. These are however unlikely to be accessible at 
the LHC, especially in the rare 6677 decay mode, although they are expected to be relevant 
at future colliders. We also explored the possibility of using jet substructure techniques to 
improve our analysis strategy. We found that these methods can be relevant in the high 
invariant-mass tail of the rrihh, distribution, namely for rrihh > 1 TeV. The impact of such 
improved analysis seems to be marginal at the LHC, whereas at future higher-energy collid¬ 


ers it can lead to a significant increase in the sensitivity on O g (see Fig. 13) and, possibly, 
on dimension -8 operators. 

There are a few directions in which our analysis could be improved, which we leave for 
future investigation. In the present work, for example, we did not try to optimize the mass 
windows for the reconstruction of the bb and 77 pairs. The choice of these windows has 
a strong impact on the non-resonant backgrounds, as they scale roughly linearly with the 
window size, and it is not unreasonable to expect that some improvement could be achieved 
with a more sophisticated analysis. Another point that we did not fully explore is the 
estimate of the backgrounds with fake photons. The analogy with single-Higgs production 
suggests that it should be possible to reduce these backgrounds to a subdominant level, but 
a thorough quantitative analysis is required to clarify this issue. 

As previously remarked, in this paper we were mostly interested in showing how double 
Higgs searches can be performed and interpreted in the general framework of effective held 
theory. We thus adopted a very simple analysis strategy avoiding unnecessary complications. 
It is clear, however, that the use of a more advanced procedure, as for instance a multivariate 
analysis, could be useful to improve the final sensitivity. In view of future high-energy collid¬ 
ers, it would also be interesting to investigate how to extract information about dimension -8 
operators. An possible strategy, in this case, could be the use of angular distributions and 
a more systematic analysis of the high-energy tails of kinematic distributions. 

Finally, a major point that is still missing in the literature is a full analysis combining 
the results obtained in the various double Higgs channels. Although 77 bb is undoubtedly 
the easiest final state to look for, other channels with larger rates can be competitive and 
might lead to a higher sensitivity on the BSM coefficients. The estimates vary significantly 
in the literature depending on the details of the simulation and of the statistical analysis. 
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The highest sensitivity is expected to come from the bbr + r~ final state, which has been 
claimed to lead to a determination of C 3 at the 40 — 60% level at the high-luminosity LHC 
with a/s = 14TeV and L = 3ab _1 [5J El [131 154] . A slightly worse sensitivity on c 3 has been 
estimated for the bbWW* channel [SII54]. while extracting the trilinear coupling from the bbbb 
channel will be much harder. For the latter case it has been estimated that the signal can 
be distinguished from a background-only hypothesis with 95% probability if c 3 < 1.2 m, 
as smaller values of c 3 lead to larger cross sections. For a fair comparison with our results, 
it should be noted that none of the above estimates, with the exception of those of Ref. [13J, 
includes the uncertainty implied by the presence of other effective operators besides the one 
controlling the trilinear coupling. Our analysis suggests that this will have an important 
impact in all cases in which a high precision can be reached on the trilinear coupling. 

All analyses of double Higgs production, including the one illustrated in this work, agree 
in concluding that the observation of double Higgs production will be a challenge for the 
LHC, even at its high-luminosity upgrade. Extracting the Higgs trilinear coupling with 
sufficient precision will thus most certainly require combining as many different channels as 
possible. This will also open up the possibility of testing precisely New Physics by following 
the strategy outlined in this paper. 

Note added: During the completion of this work Ref. appeared, where an analysis of 
gg —» hh —» 6677 is performed at a future 100 TeV collider in the context of the Standard 
Model. The authors estimate a ~ 40% sensitivity on the Higgs trilinear coupling with 3 ab -1 
of integrated luminosity, which roughly agrees with our result. Their simulation differs in 
several aspects compared to ours: on the one hand backgrounds with fake photons are 
included and their importance is pointed out; on the other hand backgrounds are computed 
at LO (including the 6677 continuum, whose simulation does not include extra jets at the 
ME level) and jet substructure techniques are not used. 
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Appendix A: Angular decomposition 


By angular momentum conservation, the scattering amplitude can be decomposed into 
two contributions, M 0 and M 2 , mediating respectively J z = 0 and J z = ±2 transitions [158]: 

A(g(Pa)g(Pb) -A h(p c )h(p d )) = {Pq V Mq + P^M 2 ) e^{p a )e v {p b ). 
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and py = 2 (p a Pc)(PbPc) / (PaPb) ~p\ is the transverse momentum of the Higgses. The expres¬ 
sion of Mq and M 2 is given by 
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where the form factors F a , Fa and G n are given in Ref. [59], and c 9 do, c 9 d 2 are the coefficients 
of the higher-derivative operators 


9* 


87 r 2 u 2 m^ . 


c 9 do + c sD2 {rrdptf&h -4»)G^(? 


a ct 


(A4) 


which give a next-to-leading correction to the non-linear Lagrangian Q. They are related 
to the coefficients of the dimension -8 operators in ([7]) by the following simple relations: 
CgDO = (47r/a 2 )c 9 2XI, CgD2 = (47 T/a 2 )c g D2- 


Appendix B: Simulation 

1. Cross sections of signal and backgrounds 

The signal and background cross sections at the generation level are summarized in 
Table |xj We simulate one extra jet from the matrix-element in addition to 6677 system in 
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hh ( 7766 ) 

77 bb 

77 33 

tth (W 77 ) 

bbh ( 6677 ) 

zh ( 6677 ) 

14 TeV 

9.71 xlO ” 2 

192 

40.3 xlO 2 

1.42 

0.14 

0.23 

100 TeV 

4.04 fb 

579 

10.5 xlO 3 

86.4 

1.47 

2.67 


TABLE X: Cross sections (in units of fir) for signal and backgrounds at 14TeV and 100 TeV after 
generation-level cuts (see text for details). 


all backgrounds except tth and 77 jj. No generation-level cuts have been applied to the 
resonant backgrounds tth and Zh. The bbh samples are restricted to satisfy prib) > 20 GeV 
at 14 TeV (increased to 30 GeV at 100 TeV) within \rj(b) | < 3 as well as A R(b,b) > 0.3. 
All extra partons are required to be within 1 77 1 <5 and their p^ threshold are set to xqcut 
in MadGraph5. The size of 77 jj is reduced by demanding pr(j) > 15 GeV at 14 TeV 
(25 GeV at 100 TeV) within \rj(j)\ < 5 and pr( 7 ) > 20 GeV (30 GeV at 100 TeV) within 
( 77 ( 7)1 < 2.5. Jets and photons are required to be separated in A R by 0.3. Appropriate 
branching fractions of Z —> bb (15.12%) and h —1 77 (0.228%), h —* bb (57.7%) are folded 
in the cross sections in Table |Xl 


2. NLO estimate of 77 bb background and k-factor 


As mentioned in Section IV, we found that the cross section of the 77 bb background 


is significantly enhanced (roughly by a factor 2 ) after including diagrams with one extra 
parton at the matrix-element level. Considering that this is the most important background 
included in our analysis it is worth exploring more in detail the origin of such enhancement. 
No analytic NLO calculation is available in the literature, so any study should currently 
rely on the use of numerical Montecarlo computations. In addition to the matched sample 
used for our analysis, we generated two additional ones. The first is obtained through a 
full NLO simulation performed with MADGRAPH5_aMC@NLO v2.1.1, including both real 
emissions and virtual corrections; it has been showered through PYTHIA 8 . The second 
is a LO simulation without extra partons at the ME level, showered through PYTHIA 6 
(with the choice of ^-ordered shower scheme). The three sets of samples are compared in 
Table |Xl| While the same version of MADGRAPH5_aMC@NLO has been used to generate all 
three samples, different sets of generation cuts had to be imposed (partly as a consequence 
of different matching procedures), making the comparison subtle. In fact, cuts have been 
imposed on “jets” (defined to be the clustered partons before parton shower) in the case 
of the NLO sample, while they apply to partons for the LO ones. For the NLO sample we 
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y/s = 14 TeV 

7766 at LO, matched up to rij jet 

7766 at NLO 

nj = 0 

nj = 1 

a [fb] at generation 

111 fb 

192 fb 

287 (223) fb 

Cuts in Eq. ( 

23 

) 

3619 

6919 

6581 

Cuts in Eq. ( 

24 

) 

557 

1274 

1164 

Cuts on m \| co 


148 

334 

301 

Cuts on m 7 e 7 ° 


13.7 

24.2 

24.5 


TABLE XI: Cut flow of the three 6677 samples described in the text at 14TeV. The number in 
parentheses corresponds to the NLO cross section when the parton-level events are restricted to 
the same mass window as in the last line of Eq. (B2). 


e 7 = 1 


imposed (see Eq. (3.5) of Ref. [38] for a definition of e 7 and n) 28 , 

Ranti-kt = 0.3 , p T (j) > 15 GeV , \v(j)\ < 5 , 

p T (b)> 20 GeV , |t 7 ( 6 )| < 3 , 

Pt(t) > 20 GeV , ( 77 ( 7 )! < 2.5 , R iS o = 0.3 . 

whereas the following cuts (with QCUT = 7 GeV) 

Pr(j) > 3 GeV , \r){j)\ < 5 , 

p T (b)> 20 GeV, ( 77 ( 6)1 < 3 , 

AR(b, 7 ) > 0.3 , AR(b, 6 ), Ai?( 7 , 7 ) > 0.3 , 

Pt{i) > 20 GeV , ( 77 ( 7 )! < 2.5 , R iso = 0.3 , 

30 < m bb < 300 GeV , 30 < m 77 < 300 GeV, 


(Bl) 


n = 1, 


(B2) 


e 7 = 1 


n — 1 , 


have been imposed on the LO samples. The prescription of Ref. [38] for the photon isolation 
is important especially for the NLO computation to guarantee the correct cancellation of 
IR divergences. We apply the same photon isolation to the tree-level matched samples 
for consistency 29 (choosing this option automatically inactivates the cut on A R(j, 7 ) in 
MADGRAPH5_aMC@NLO). The low value of px(j) in Eq. (B2) is due to the low value of 
the matching scale. The lower bound of the m bb mass window in Eq. (B2) has been chosen 
so as to obtain an optimal matching scale and, at the same time, sufficiently small to avoid 
possible distortions in the signal region. 


The value of the cross section for each of the three samples is reported in Table XI 


The entries in the first row should not be literally compared, since they refer to samples 


28 We thank Marco Zaro for helping us to impose appropriate cuts on 6-jets in our NLO simulation. 

29 However, we find no practical discrepancy between the step-function type cone isolation, i.e. AR("/,X) > 
0.3 where X denotes any object near the photon, and the prescription in |38j in the tree-level matched 
samples. 
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generated with different sets of cuts and procedures. These generation-level differences, 
however, become irrelevant after imposing the first round of cuts leading to the number of 


events in the second row of Table XI The k-factor, i.e. the ratio of the fourth to second 
column, is approximately 2 in each step of cut-flow. A similar k-factor is also found by 
taking the ratio of the third to second column, i.e. of LO samples with and without one 
extra parton. This shows that, at least at the level of cross section, the matched sample gives 
a very good approximation of the NLO one. This indicates that virtual corrections give a 
small contribution and the bulk of the NLO correction comes from real emissions. Another 
important aspect are the differential distributions. Although total rates agree well, the shape 
of distributions could differ between the NLO sample and the LO matched one, leading to 
different cut efficiencies. Figure [2T| compares the distributions of several kinematic variables 
relevant for our analysis. While the higher value bins are subject to statistical fluctuations, 
overall the differential shapes of the two samples agree well. This justifies our use of the 
LO matched 77 bb sample -whose generation requires much less CPU time- throughout our 
analysis. We expect that the same results will also apply in the case of a 100 TeV collider. 
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